Activity in Pectoral Sandpipers: Statistics and figures
Background
Purpose: This notebook is used to perform statistical analyses and create tables/figures that are included in the manuscript.
Inputs: The following input datasets are required:
- Capture data (data/capture-data.csv)
- GPS data (data/gps-data.csv)
- Nest data (data/nest-data.csv)
- Parentage data (data/parentage-data.csv)
- Sighting data (data/sighting-data.csv)
- Cleaned ODBA data with activity thresholds (outputs/processed-data/ODBA-and-activity_10min-intervals_onboard.csv)
- Male trips offsite (outputs/processed-data/male-trips-offsite.csv)
Data processing steps are outlined in the next section.
Outputs: Tables and figures for the manuscript are saved in “outputs/figures-and-tables”
Notes on data processing steps
Prior to performing the analyses below, steps for cleaning and preparing the data are performed in separate scripts:
- clean-onboard-odba-data.R: Cleans data processed on board the Druid accelerometer.
- define-activity-thresholds.Rmd: Determines activity thresholds for each individual
- map-male-locations.Rmd: Identifies when males took trips outside the study area
These processing scripts are found in the “scripts” folder. The key outputs are stored in “outputs/processed-data”.
Import data and prepare workspace
Most packages used in this Notebook can be installed using the install.packages function. One exception is “wadeR” which was developed by researchers at MPIO and needs to be installed from github.
# library(remotes)
# remotes::install_github("mpio-be/wadeR") # Utqiagvik (Barrow) mappinglibrary(tidyr) # for reshaping datasets
library(dplyr) # for other data manipulation
library(ggplot2) # for plotting
library(scales) # for plotting dates
library(lubridate) # for working with dates and times
library(ggalt) # for dumbbell charts
library(data.table) # for fast table reading and rolling functions
library(flextable) # generates neat tables
library(officer) # formatting flextable for Microsoft Word
library(fasttime) # for faster conversions to posixct
library(wadeR)
library(stringr) # count matches in string
library(here) # for easier directory navigation
library(glue) # for easier naming and referring to directories
library(dichromat) # for checking colours in figures
library(emmeans) # for computing estimated marginal means (EMMs)
library(effects) # for visualising effects
library(glmmTMB) # for glmm with beta error distribution
library(geosphere) # for distances lat lon
library(lme4) # for mixed effects models
library(spatialrisk) # for distances lat lon
library(patchwork) # better than ggarrange for arranging figures
library(DHARMa) # for checking model assumptions
library(sf) # for geometry functions
library(broom.mixed) # for converting glmmtmb to flextableSome functions are used in multiple scripts and are therefore stored in their own script for easier access (general-helper-functions.R). This helper script needs to be loaded for the analyses below to run.
source("./scripts/general-helper-functions.R")Male activity was recorded using miniature accelerometers and received in two different ways: processed on board the device (“onboard” data) or stored in raw form (“rawbased” data). During the course of the study, we also considered different methods for processing the data. These included summarising the data at either 1-min or 10-min intervals, and either defining activity thresholds for each individual separately (“individual” thresholds) or together (“global” thresholds).
There are advantages and disadvantages of different approaches, but results presented in the main text of the manuscript are based on the “onboard” data with a window-length of 10 minutes and “individual” thresholds.
In the supplementary material of the manuscript, we present a comparison of these different approaches. By changing the definitions in the code chunk below, and running the appropriate pre-processing steps, analyses in this R Notebook can be performed using datasets processed in different ways. These alternative datasets and pre-processing scripts are available upon request.
window_length = 10 # 1 or 10
data_type = "onboard" # rawbased or onboard
threshold_type = "individual" # global or individual
if(!(window_length==1 | window_length==10)){
print("Error: Select a window length of either 1 or 10")
}
if(!(data_type=="onboard" | data_type=="rawbased")){
print("Error: Select a data_type of either 'onboard' or 'rawbased'")
}
if(!(threshold_type=="individual" | threshold_type=="global")){
print("Error: Select a threshold type of either 'individual' or 'global'")
}
if(threshold_type=="global" & window_length!=10){
print("Warning: We recommend using a window length of 10 when using the global threshold. The global threshold was developed using data summarised at 10-minute intervals.")
}
if(threshold_type=="global" & data_type!="onboard"){
print("Warning: We recommend using the 'onboard' data when applying the global threshold. The global threshold was developed using data processed on board the device.")
}
if(window_length=="1" & data_type!="onboard"){
print("Error: Onboard data is only available with a window length of 10 minutes.")
}
if(threshold_type == "global") {
global_threshold = 0.27 # this threshold is derived from a separate analysis that is available upon request
}# 10-min onboard activity data is used for the main analyses in the manuscript
# all other data types are stored in a supplementary folder
if(window_length==10 & data_type=="onboard") {
activity_data_path = here("outputs","processed-data")
} else {
activity_data_path = here("supplementary", "supp-data")
}captures <- read.csv(here("data","capture-data.csv"))
nests <- read.csv(here("data","nest-data.csv"))
sightings <- read.csv(here("data", "sighting-data.csv"))
parentage <- read.csv(here("data", "parentage-data.csv"))
gps <- read.csv(here("data","gps-data.csv"))
activity_onboard <- read.csv(here("outputs", "processed-data",
"ODBA-and-activity_10min-intervals_onboard.csv"))
# this dataset is used to help determine when final data transmission occurred
activity <- read.csv(glue("{activity_data_path}/ODBA-and-activity_{window_length}min-intervals_{data_type}.csv"))
# this dataset is the source of ODBA and activity data for analyses
# there is the option to customise the input to analyse different activity datasets
male_trips_offsite <- read.csv(here("outputs", "processed-data", "male-trips-offsite.csv"))All dates and times are local Alaskan times, unless specified otherwise.
captures <- captures %>%
mutate(date=local_date_to_posix(as.character(date)),
released_dt = local_dt_to_posix(paste(date,released)))
nests <- nests %>%
mutate(date_found=local_date_to_posix(as.character(date_found)),
est_hatch_date=local_date_to_posix(as.character(est_hatch_date)),
hatch_date=local_date_to_posix(as.character(hatch_date)))
sightings <- sightings %>%
mutate(datetime=local_dt_to_posix(as.character(datetime)))
gps <- gps %>%
mutate(datetime=local_dt_to_posix(as.character(datetime)))
activity_onboard <- activity_onboard %>%
mutate(window_end_utc=acc_to_posix(as.character(window_end_utc)),
window_end_localtime=dt_to_local(window_end_utc))
activity <- activity %>%
mutate(window_end_utc=acc_to_posix(as.character(window_end_utc)),
window_end_localtime=dt_to_local(window_end_utc))
male_trips_offsite <- male_trips_offsite %>%
mutate(start_dt_local=local_dt_to_posix(as.character(start_dt_local)),
end_dt_local=local_dt_to_posix(as.character(end_dt_local)))To avoid biases that might arise from incomplete data, we use a minimum threshold for data availability for some analyses. Specifically, we have defined a minimum number of hours of data that need to be available per individual per day (“min_acc_hrs_per_day”) and a minimum number of days of data that need to be available per individual (“min_acc_days”). Data that do not meet these requirements (i.e. insufficient data available from an individual on a particular day, or insufficient data available from a particular individual overall) are excluded from some analyses. By default, we have required at least 85% complete data per individual per day, and at least 2 days of data from each individual.
min_acc_hrs_per_day = 20.4 # 20.4 means at least 85% complete
min_acc_days = 2 # except where otherwise specifiedmale_colour = "darkgreen"
female_colour = "gray32"
success_colour = "#2267E7" #"#20EAABFF"
fail_colour = "#DF536B"
default_shape = 21
uncertain_shape = 4
legend.heading.size = 10
axis.heading.size = 10
label.size = 10
font.family = "Calibri"
small.point.size = 1
# ribbon plot specifications
point_size = small.point.size+1
point_alpha = 0.5
col_scale <- c("springgreen4", "darkgray")
extra_nest_shape = 13
windowsFonts("Calibri" = windowsFont("Calibri")) # avoids font issuesset_flextable_defaults(
font.size = 11,
theme_fun = theme_booktabs,
font.family = "Calibri",
padding = 3,
background.color = "white")
sect_properties <- prop_section(
page_size = page_size(orient = "portrait",
width = 21 * 0.393701,
height =29.7 * 0.393701),
type = "continuous",
page_margins = page_mar(bottom = 2.54 * 0.393701,
top = 2.54 * 0.393701,
left = 2.54 * 0.393701,
right = 2.54 * 0.393701)
)save_plots = TRUE
save_tables = TRUE
output_path <- here('outputs','figures-and-tables')
output_supp_effects <- here('supplementary','supp-outputs','effect-sizes')
if (!dir.exists(output_path)) {
dir.create(output_path)
}Preparing data for analysis
Note: The tables and figures in this section are for context and data-checking only. They are not included in the manuscript.
Identify tagged males
tagged_males <- captures %>%
filter(tag_action=="on") %>%
filter(!colour_combo %in% colour_combos_to_exclude) %>%
select(colour_combo, bird_ID, tag_ID, date, released_dt, body_mass)ggplot(tagged_males, aes(x=date))+
theme_classic()+
geom_bar(stat="count",
fill=male_colour) +
scale_x_datetime(breaks=date_breaks(width="1 day"),
date_labels="%d",
expand = c(0, 0)) +
coord_cartesian(xlim=c(local_date_to_posix("2022-06-06"),
local_date_to_posix("2022-06-28")))+
xlab("\n Date (June)")+
ylab("Number of new tags deployed on males\n")
In 2022, we equipped 100 adult male pectoral sandpipers with 3.7 g accelerometer tags (weight includes harness and glue) with GPS function (Debut NANO, Druid Inc.) between 2022-06-07 and 2022-06-20. We targeted males who were exhibiting territorial behaviours. We stopped tagging males when we ran out of tags (although new males were still entering the study site).
The weight of the device was between approximately 3.2 and 4.2 % of the bird’s body mass.
Determing female fertility
We defined the fertile period of each female pectoral sandpiper as the interval beginning 5 days before the female began laying eggs (Kempenaers and Valcu 2017 Nature, https://doi.org/10.1038/nature20813) and ending the day the penultimate egg was laid. Clutch size ranged from 2 to 4 eggs (89.7% of nests contained 4 eggs, 6.9% contained 3 eggs, and 3.4% contained 2 eggs), meaning that the fertile period for each female ranged from 6-8 days. It is not actually known whether females still engage in copulations after initiating their clutch.
If the nest was found during laying (clutch size when found < maximum clutch size), the beginning of the fertile period was defined as date found - clutch size when found (females lay one egg per day).
If the nest was complete when found, the beginning of the fertile period was defined as: hatch date - the typical incubation period (22 days) - (maximum clutch size - 1) . For hatch date, we used either the earliest actual hatch date if the eggs successfully hatched in the incubator, or the earliest predicted hatch date (based on the floatation method) if they did not.
Our calculations assume that whenever the nest was found or checked, the female had not laid any eggs yet on that day.
nests$found_status <- if_else(
nests$first_clutchsize==nests$max_clutchsize,
"complete",
"laying"
)
nests <- nests %>%
mutate(clutch_init = case_when(found_status=="laying" ~ date_found-first_clutchsize*(24*60*60),
is.na(hatch_date) ~ est_hatch_date - 22*(24*60*60) - (max_clutchsize-1)*(24*60*60),
TRUE ~ hatch_date - 22*(24*60*60) - (max_clutchsize-1)*(24*60*60)),
penult_egg = case_when(max_clutchsize>1 ~ clutch_init + (max_clutchsize-2)*(24*60*60),
TRUE ~ clutch_init),
first_fertile_day=clutch_init-5*(24*60*60),
last_fertile_day = penult_egg)For the overall study site, we defined the start of the fertile period as the 5% quantile of all fertile start dates, and the end as the 95% quantile. A similar approach was used to define the end of the fertile period in Lesku et al. 2012 Science (https://doi.org/10.1126/science.1220939.
After the start date is ‘fertile’, end date onwards is ‘post-fertile’.
start_FertilePeriod <- floor_date(
quantile(nests$first_fertile_day, 0.05, na.rm=T),
unit="day")
end_FertilePeriod <- floor_date(
quantile(nests$last_fertile_day, 0.95, na.rm=T),
unit="day")all_FertileDates <- seq(min(nests$first_fertile_day, na.rm=T),
by = "day",
length.out = max(nests$last_fertile_day, na.rm=T)-
min(nests$first_fertile_day, na.rm=T))
countFertileFemales <- function(date) {
count <- date>=nests$first_fertile_day &
date<=nests$last_fertile_day
sum <- sum(count, na.rm=TRUE)
return(sum)
}
ff_by_date <- data.frame(date=unlist(all_FertileDates))
ff_by_date$n_FertileFemales <- sapply(all_FertileDates,
countFertileFemales)
barplot_ff <- ggplot(ff_by_date, aes(x=date,
y=n_FertileFemales)) +
# theme_classic()+
theme(panel.background=element_blank(),
axis.line = element_line(size=0.5, colour="black"))+
annotate(geom="rect",
xmin=start_FertilePeriod-24*60*30,
xmax=end_FertilePeriod+24*60*30,
ymin=0,
ymax=43,
fill="grey",
alpha=0.3)+
geom_bar(stat="identity",
fill=female_colour) +
coord_cartesian(xlim=c(local_date_to_posix("2022-06-06"),
local_date_to_posix("2022-07-02")))+
scale_y_continuous(limits=c(0,43), expand=c(0,0))+
xlab("\n Date") +
ylab("Number of fertile females\n")
barplot_ff
Summarise paternity data
We collected blood samples from all of the males we tagged (as well as some that we didn’t) and from eggs/chicks on our study plot. DNA was extracted from these samples by collaborators at the University of Alaska. These DNA samples were shipped to the research institute.
parentage <- left_join(parentage, captures[,c("bird_ID","colour_combo")],
by=join_by("ID.father"=="bird_ID")) %>%
mutate(ID.father=as.character(ID.father)) %>%
mutate(ID.father=replace_na(ID.father,"unbanded")) %>%
rename(LL.father=colour_combo) %>%
mutate(LL.father=replace_na(LL.father,"unbanded")) %>%
left_join(., captures[,c("bird_ID","colour_combo")], by=join_by("ID.mother"=="bird_ID")) %>%
rename(LL.mother=colour_combo)fathers <- parentage %>%
group_by(ID.father, LL.father, Nest) %>%
summarise(n_eggs=n()) %>%
group_by(ID.father, LL.father) %>%
summarise(n_eggs=sum(n_eggs),
n_nests=n())
fathers %>% flextable()ID.father | LL.father | n_eggs | n_nests |
|---|---|---|---|
141269009 | PI,PI,R | 6 | 2 |
141269026 | O,DB,O | 4 | 1 |
141269043 | DG,O,O | 4 | 1 |
141269044 | Y,DG,PI | 5 | 2 |
141269046 | DB,DG,Y | 4 | 1 |
141269114 | PI,R,PI | 3 | 1 |
141269125 | Y,Y,R | 1 | 1 |
141269127 | DG,DB,DG | 4 | 1 |
141269144 | DB,Y,R | 3 | 1 |
141269146 | R,DG,PI | 2 | 1 |
141269147 | O,DG,DG | 5 | 2 |
141269148 | DB,R,O | 3 | 1 |
141269152 | Y,Y,Y | 1 | 1 |
141269169 | DB,R,R | 4 | 1 |
141269175 | DG,R,R | 3 | 1 |
141269182 | O,Y,O | 2 | 1 |
unbanded | unbanded | 136 | 38 |
(Table for background only; not to be included in paper)
parentage %>%
group_by(Nest) %>%
summarise(fathers=toString(unique(LL.father)),
unique_fathers=n_distinct(LL.father)) %>%
filter(fathers!="NA" & fathers!="unbanded") %>%
arrange(desc(unique_fathers)) %>%
flextable()Nest | fathers | unique_fathers |
|---|---|---|
P1103 | unbanded, DG,R,R | 2 |
P1401 | unbanded, Y,Y,Y | 2 |
P201 | unbanded, Y,DG,PI | 2 |
P208 | unbanded, O,Y,O | 2 |
P403 | PI,R,PI, O,DG,DG | 2 |
P701 | R,DG,PI, PI,PI,R | 2 |
P702 | Y,Y,R, DB,Y,R | 2 |
P1001 | DG,DB,DG | 1 |
P1003 | O,DB,O | 1 |
P1403 | Y,DG,PI | 1 |
P1903 | DB,DG,Y | 1 |
P1905 | DG,O,O | 1 |
P203 | DB,R,O | 1 |
P204 | DB,R,R | 1 |
P214 | O,DG,DG | 1 |
P401 | PI,PI,R | 1 |
Exclude unreliable or irrelevant GPS data
GPS data are filtered down to data collected while the device was on the bird.
Occasionally, the devices recorded a latitude and longitude of 200. This is an error value that doesn’t mean anything and should be excluded.
male_gps <- gps %>%
left_join(., tagged_males, join_by('tag_ID')) %>%
filter(lat!=200 & lon!=200) %>% # these are nonsense values
filter(datetime>released_dt) %>%
mutate(local_date=floor_date(datetime, unit="days"),
bird_ID=as.factor(bird_ID)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE) Determining site tenure
We use capture date as a proxy for arrival date.
We use the latest of sighting, GPS and ODBA data as a proxy for departure.
There are two birds who were “seen” after they last transmitted data. One of these was definitely a mistake; the wrong GPS point must have been recorded, because the sighting is associated with a video file from 10 days earlier (14 June rather than 24 June). So for this bird (DG,DB,DG) I am using the last date of data transmission as his departure date.
The other I have no reason to exclude, and the sighting was only 1 day after he last transmitted data.
latest_onboard = activity_onboard %>%
group_by(tag_ID) %>%
mutate(window_end_localtime_hour=floor_date(window_end_localtime, unit="hours")) %>%
summarise(last_ODBA = max(window_end_localtime, na.rm=T),
last_hour_n_ODBA = sum(window_end_localtime>=(last_ODBA-60*60), na.rm=T))
latest_gps = male_gps %>%
group_by(tag_ID) %>%
summarise(last_GPS = max(datetime, na.rm=T),
last_hour_n_GPS = sum(datetime>=(last_GPS-60*60), na.rm=T))
latest_sightings <- inner_join(sightings, tagged_males, by='colour_combo') %>%
group_by(tag_ID) %>%
summarise(last_sighting = max(datetime, na.rm=T),
last_hour_n_sightings = sum(datetime>=(last_sighting-60*60), na.rm=T))
male_departures = left_join(tagged_males, latest_onboard, join_by("tag_ID")) %>%
left_join(., latest_gps, join_by("tag_ID")) %>%
left_join(., latest_sightings, join_by("tag_ID"))
male_departures = male_departures %>%
mutate(last_data_akdt = pmax(last_ODBA,last_GPS,last_sighting, na.rm=T),
last_data_utc = dt_to_acc(last_data_akdt),
last_data_cet = with_tz(last_data_akdt, tz="CET"),
across(c(last_ODBA, last_GPS, last_sighting),
~case_when(. == last_data_akdt ~ cur_column()), .names = 'new_{col}')) %>%
unite(col=last_data_type, starts_with('new'), na.rm = TRUE, sep = ', ') %>%
mutate(last_hour_n = case_when(last_data_type=="last_ODBA" ~ last_hour_n_ODBA,
last_data_type=="last_GPS" ~ last_hour_n_GPS,
last_data_type=="last_sighting" ~ last_hour_n_sightings,
TRUE ~ NA))
male_tenure = male_departures %>%
select(colour_combo, bird_ID, tag_ID, released_dt, last_data_akdt, last_data_type, last_hour_n) %>%
rename(arrival_dt=released_dt,
departure_dt=last_data_akdt) %>%
mutate(arrival_date=floor_date(arrival_dt, unit="days"),
departure_date=floor_date(departure_dt, unit="days"),
departure_dt_confidence = case_when(last_hour_n>1 & last_data_type!="last_sighting" ~ "OK",
TRUE ~ "low"),
unrevised_tenure_days_exact = as.numeric(difftime(departure_dt, arrival_dt, unit="days")),
unrevised_tenure_days = as.numeric(difftime(departure_date, arrival_date, unit="days")))
rm(male_departures, latest_onboard, latest_gps, latest_sightings)Some males left and then returned again. We subtract the days that they were away from tenure.
male_trips_offsite <- male_trips_offsite %>%
mutate(offsite_days = as.numeric(difftime(as.Date(end_dt_local),
as.Date(start_dt_local),
unit="days")))
male_trips_offsite %>%
summarise(mean_offsite_days=mean(offsite_days),
sd_offsite_days=sd(offsite_days),
min_offsite_days=min(offsite_days),
max_offsite_days=max(offsite_days))## mean_offsite_days sd_offsite_days min_offsite_days max_offsite_days
## 1 3.57 2.07 1 7male_offsite_summary = male_trips_offsite %>%
select(bird_ID, duration_hrs, offsite_days) %>%
rename(offsite_hrs=duration_hrs)male_tenure <- left_join(male_tenure, male_offsite_summary, join_by('bird_ID')) %>%
mutate(tenure_days_exact=case_when(is.na(offsite_hrs) ~ unrevised_tenure_days_exact,
TRUE ~ unrevised_tenure_days_exact-(offsite_hrs/24)),
tenure_days=case_when(is.na(offsite_days)~unrevised_tenure_days,
TRUE ~ unrevised_tenure_days-offsite_days)) %>%
select(colour_combo, bird_ID, tag_ID,
arrival_dt, departure_dt, arrival_date, departure_date,
tenure_days, tenure_days_exact)Examine breeding density relative to other years
In total, there were 58 breeding females and median male tenure was 10 days. This was therefore a moderate density breeding season, compared to other years (Kempenaers and Valcu 2017, Nature; number of breeding females ranges from approx. 23-90). It was higher density year than 2008 and 2009, which were the breeding seasons when data were collected for the previous activity study (Lesku et al. 2012, Science: 39 nests in 2008 and 17 nests in 2009.)
Exclude unreliable or irrelevant ODBA data
A simple and common measure of activity is overall dynamic body acceleration (ODBA). ODBA is a measure of total movement, combined across the three axes of the accelerometer. It is calculated by first subtracting static acceleration from each axis (acceleration related to posture; most often calculated by taking the 2-second running mean) to get dynamic acceleration, then summing together the absolute dynamic acceleration for each axis (for example, see: https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00057.x) ODBA is also commonly used as a proxy for energy expenditure related to movement.
Data have already been excluded when the tag was not on the bird. In a separate script ( define-activity-thresholds.RmD ) I also defined an activity threshold for each individual bird.
Below, I am excluding the following data:
- Activity thresholds from birds without reliable estimates (due to <36 hours total available data)
- Data collected when birds were offsite
- Data collected <12 hours after putting the tag on
min_hrs_after_release = 12# Exclude activity thresholds with <36 hours of data (when using individual threshold)
if(threshold_type=="individual") {
activity <- activity %>%
mutate(activity_threshold=case_when(n_windows_threshold<=36*60/window_length ~ NA,
TRUE ~ activity_threshold),
active=case_when(mean_ODBA>activity_threshold ~ "active",
mean_ODBA<=activity_threshold ~ "inactive",
TRUE ~ NA))
}
# Exclude trips offsite
activity <- male_trips_offsite %>%
select(colour_combo, bird_ID, start_dt_local, end_dt_local) %>%
right_join(., activity, join_by('colour_combo','bird_ID'),
relationship="one-to-many") %>%
filter(is.na(start_dt_local) |
window_end_localtime < start_dt_local |
window_end_localtime-window_length*60 > end_dt_local) %>%
select(!start_dt_local & !end_dt_local)
# Exclude 12 hours after release
activity <- tagged_males %>%
select(colour_combo, bird_ID, released_dt) %>%
right_join(., activity, join_by('colour_combo','bird_ID')) %>%
filter(window_end_localtime-window_length*60 > released_dt+(min_hrs_after_release*60*60))There was also one male who travelled 300 km in one day (out over the ocean for 150km, then back; see map-male-locations.Rmd). It doesn’t make sense to include this data in our analysis. He was R,Y,DB and his travel date was 20 June.
activity <- activity %>%
mutate(window_end_local_date=floor_date(window_end_localtime, unit="day"))
activity <- activity %>%
filter(!(colour_combo=="R,Y,DB" & window_end_local_date==as.Date("2022-06-20")))Summarise and define activity thresholds and levels
if(threshold_type=="individual") {
activity %>%
summarise(min_threshold=min(activity_threshold, na.rm=T),
max_threshold=max(activity_threshold, na.rm=T),
median_threshold=median(activity_threshold, na.rm=T)) %>%
flextable()
} else {
print(paste0("Global threshold = ", global_threshold))
}min_threshold | max_threshold | median_threshold |
|---|---|---|
0.147 | 0.37 | 0.239 |
By determining an activity threshold for each individual, based on the distribution of his activity data, we should minimise bias in activity estimates caused by variation in measurements by different devices and minor variations in attachment. However, it is possible that we are also sometimes underestimating the activity of highly active birds and overestimating the activity of less active birds, since these distributions in activity distributions can affect the values of the threshold. Compared to having a global threshold across all birds, this method is therefore more conservative and more likely to underestimate inter-individual variation in activity levels.
Note: Birds without activity thresholds (i.e. <36 hours data) have not been excluded from the “activity” dataset; for some analyses, the ODBA values of these birds are still relevant and worth including.
if(threshold_type=="global") {
activity <- activity %>%
mutate(activity_threshold=global_threshold,
active=case_when(mean_ODBA>global_threshold ~ "active",
TRUE ~ "inactive"))
}Calculate daily activity
activity$local_date = floor_date(activity$window_end_localtime, unit="days")
daily_activity <- activity %>%
#filter(!is.na(activity_threshold)) %>%
group_by(bird_ID, tag_ID, colour_combo, local_date) %>%
summarise(activity_threshold=first(activity_threshold),
n_windows=n(),
n_hours=n_windows*window_length/60,
perc_active=sum(active=="active")/n_windows*100,
total_active=sum(active=="active"),
mean_ODBA=mean(mean_ODBA))Determine distance of males from fertile nests
predated_nests <- nests %>%
filter(is.na(est_hatch_date)) %>% # predated nests had no estimated hatch date, eggs not floated
select(nest, date_found, first_clutchsize, lat, lon)
predated_nests %>% flextable()nest | date_found | first_clutchsize | lat | lon |
|---|---|---|---|---|
P1109 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 |
P1110 | 2022-06-30 00:00:00 | 3 | 71.3 | -157 |
P1111 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 |
P1112 | 2022-07-01 00:00:00 | 2 | 71.3 | -157 |
P1409 | 2022-06-28 00:00:00 | 4 | 71.3 | -157 |
P1910 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 |
P206 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 |
P207 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 |
To estimate the fertile dates of predated nests, I’ve assumed that the maximum clutch size would have been 4 (so nests found with clutch size < 4 were found during laying).
For nests that were not found during laying, we only know that the female was no longer fertile when the nest was found. All of the nests were found after the end of the fertile phase for the plot (based on the fertile periods of other nests). I am therefore using the entire fertile phase for the plot as the fertile period for these predated nests. This is used to help establish which males were spending time near the nest when it might have been fertile.
predated_nests <- predated_nests %>%
mutate(clutch_init = case_when(first_clutchsize<4 ~ date_found-first_clutchsize*86400,
TRUE ~ NA),
penult_egg = case_when(is.na(clutch_init) ~ NA,
TRUE ~ clutch_init + 2*86400),
first_fertile_day = case_when(!is.na(clutch_init) ~ clutch_init-5*86400,
TRUE ~ start_FertilePeriod),
last_fertile_day = case_when(!is.na(penult_egg) ~ penult_egg,
TRUE ~ end_FertilePeriod)
)
predated_nests %>% flextable()nest | date_found | first_clutchsize | lat | lon | clutch_init | penult_egg | first_fertile_day | last_fertile_day |
|---|---|---|---|---|---|---|---|---|
P1109 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 | 2022-06-10 00:00:00 | 2022-06-26 00:00:00 | ||
P1110 | 2022-06-30 00:00:00 | 3 | 71.3 | -157 | 2022-06-27 00:00:00 | 2022-06-29 00:00:00 | 2022-06-22 00:00:00 | 2022-06-29 00:00:00 |
P1111 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 | 2022-06-10 00:00:00 | 2022-06-26 00:00:00 | ||
P1112 | 2022-07-01 00:00:00 | 2 | 71.3 | -157 | 2022-06-29 00:00:00 | 2022-07-01 00:00:00 | 2022-06-24 00:00:00 | 2022-07-01 00:00:00 |
P1409 | 2022-06-28 00:00:00 | 4 | 71.3 | -157 | 2022-06-10 00:00:00 | 2022-06-26 00:00:00 | ||
P1910 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 | 2022-06-10 00:00:00 | 2022-06-26 00:00:00 | ||
P206 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 | 2022-06-10 00:00:00 | 2022-06-26 00:00:00 | ||
P207 | 2022-06-30 00:00:00 | 4 | 71.3 | -157 | 2022-06-10 00:00:00 | 2022-06-26 00:00:00 |
nests <- nests %>%
filter(!nest %in% predated_nests$nest) %>%
full_join(., predated_nests,
join_by('nest', 'date_found', 'lat','lon','first_clutchsize', 'clutch_init', 'penult_egg','first_fertile_day','last_fertile_day'))nest_parentage <- parentage %>%
filter(is.na(Comments) | Comments!="could not be genotyped") %>%
select(Nest, LL.father) %>%
unique() %>%
group_by(Nest) %>%
arrange(LL.father) %>%
summarise(father1 = first(LL.father),
father2 = case_when(last(LL.father)==father1 ~ NA,
TRUE ~ last(LL.father)),
n_fathers = n()) %>%
rename(nest = Nest)nest_dates <- nests %>%
select(nest, first_fertile_day, last_fertile_day) %>%
pivot_longer(cols = first_fertile_day:last_fertile_day,
values_to = 'date') %>%
filter(!is.na(date)) %>%
group_by(nest) %>%
complete(date = seq(min(date), max(date), by='days')) %>%
select(-name)
nest_long <- left_join(nest_dates, nests, join_by('nest'), relationship='many-to-one') %>%
left_join(., nest_parentage, join_by('nest'), relationship='many-to-one') %>%
select(nest, date, first_fertile_day, last_fertile_day, lat, lon, father1, father2) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE)male_date_gps <- male_gps %>%
group_by(colour_combo, bird_ID, tag_ID, local_date, released_dt) %>%
summarise(n_gps=n(),
centroid = st_centroid(st_union(geometry)))daily_activity_gps <- activity %>%
group_by(bird_ID, tag_ID, colour_combo, local_date) %>%
summarise(activity_threshold=first(activity_threshold),
mean_ODBA=mean(mean_ODBA),
perc_active=case_when(is.na(activity_threshold) ~ NA,
TRUE ~ sum(active=="active")/(sum(active=="active")+sum(active=="inactive"))*100),
prop_activity=case_when(is.na(activity_threshold) ~ NA,
TRUE ~ perc_active/100),
total_active_windows=sum(active=="active"),
n_windows=n(),
n_hours=n()/(60/window_length)) %>%
mutate(bird_ID=as.factor(bird_ID)) %>%
filter(n_hours>min_acc_hrs_per_day) %>%
select(bird_ID, tag_ID, colour_combo, local_date, mean_ODBA, perc_active, prop_activity, n_windows) %>%
full_join(., male_date_gps, join_by('tag_ID', 'bird_ID', 'colour_combo', 'local_date'),
relationship='many-to-many') %>%
filter(!is.na(mean_ODBA))Note that this next chunk of code takes a while to run (getting distance from each male from each nest on each day; ~ 3-5 mins).
nest_male_date <- inner_join(nest_long,
daily_activity_gps,
join_by('date'=='local_date'),
relationship='many-to-many') %>%
filter(!is.na(lat) & !is.na(mean_ODBA))
nest_male_date <- nest_male_date %>%
rowwise() %>%
mutate(dist_m = st_distance(geometry, centroid))nest_male_date <- nest_male_date %>%
filter(!is.na(n_gps)) %>%
mutate(success=case_when(colour_combo==father1 ~ 1,
colour_combo==father2 ~ 1,
TRUE ~ 0),
bird_ID = as.factor(bird_ID),
date = as.factor(date),
nest = as.factor(nest))nest_male_date %>%
filter(success==1) %>%
select(nest, father1, father2, dist_m, colour_combo, success) %>%
group_by(nest, father1, father2, colour_combo, success) %>%
summarise(mean_dist_m=mean(dist_m),
min_dist_m=min(dist_m),
median_dist_m=median(dist_m),
max_dist_m=max(dist_m),
n=n()) %>%
arrange(median_dist_m) %>%
select(!c(geometry)) %>%
flextable()nest | father1 | father2 | colour_combo | success | mean_dist_m | min_dist_m | median_dist_m | max_dist_m | n | geometry |
|---|---|---|---|---|---|---|---|---|---|---|
P208 | O,Y,O | unbanded | O,Y,O | 1 | 34.6 [m] | 8.7 [m] | 33.3 [m] | 52.7 [m] | 7 | [[XY]] |
P203 | DB,R,O | DB,R,O | 1 | 38.8 [m] | 19.9 [m] | 36.0 [m] | 73.6 [m] | 7 | [[XY]] | |
P214 | O,DG,DG | O,DG,DG | 1 | 273.9 [m] | 31.5 [m] | 57.7 [m] | 1,035.8 [m] | 8 | [[XY]] | |
P1403 | Y,DG,PI | Y,DG,PI | 1 | 62.1 [m] | 62.1 [m] | 62.1 [m] | 62.1 [m] | 1 | [[XY]] | |
P1903 | DB,DG,Y | DB,DG,Y | 1 | 84.2 [m] | 47.8 [m] | 63.0 [m] | 163.1 [m] | 4 | [[XY]] | |
P1001 | DG,DB,DG | DG,DB,DG | 1 | 66.5 [m] | 66.5 [m] | 66.5 [m] | 66.5 [m] | 1 | [[XY]] | |
P403 | O,DG,DG | PI,R,PI | O,DG,DG | 1 | 69.9 [m] | 63.8 [m] | 69.5 [m] | 76.3 [m] | 3 | [[XY]] |
P1003 | O,DB,O | O,DB,O | 1 | 91.4 [m] | 79.4 [m] | 92.1 [m] | 102.8 [m] | 7 | [[XY]] | |
P204 | DB,R,R | DB,R,R | 1 | 112.5 [m] | 105.5 [m] | 111.1 [m] | 120.8 [m] | 3 | [[XY]] | |
P201 | Y,DG,PI | unbanded | Y,DG,PI | 1 | 289.1 [m] | 43.2 [m] | 118.8 [m] | 799.9 [m] | 5 | [[XY]] |
P1905 | DG,O,O | DG,O,O | 1 | 165.3 [m] | 112.1 [m] | 119.1 [m] | 383.9 [m] | 6 | [[XY]] | |
P401 | PI,PI,R | PI,PI,R | 1 | 172.5 [m] | 142.0 [m] | 162.8 [m] | 214.6 [m] | 7 | [[XY]] | |
P702 | DB,Y,R | Y,Y,R | Y,Y,R | 1 | 173.7 [m] | 148.6 [m] | 173.7 [m] | 198.7 [m] | 2 | [[XY]] |
P1401 | Y,Y,Y | unbanded | Y,Y,Y | 1 | 348.1 [m] | 167.0 [m] | 187.8 [m] | 1,350.2 [m] | 8 | [[XY]] |
P403 | O,DG,DG | PI,R,PI | PI,R,PI | 1 | 191.0 [m] | 135.2 [m] | 208.0 [m] | 222.9 [m] | 7 | [[XY]] |
P701 | PI,PI,R | R,DG,PI | PI,PI,R | 1 | 763.5 [m] | 741.1 [m] | 768.0 [m] | 777.1 [m] | 4 | [[XY]] |
P1103 | DG,R,R | unbanded | DG,R,R | 1 | 800.7 [m] | 800.7 [m] | 800.7 [m] | 800.7 [m] | 1 | [[XY]] |
P701 | PI,PI,R | R,DG,PI | R,DG,PI | 1 | 896.3 [m] | 896.3 [m] | 896.3 [m] | 896.3 [m] | 1 | [[XY]] |
In most cases (14/18), the median daily distance of males from the nest they sired was <200 metres. There was one male whose median daily distance was slightly further (208 m), but he spent a large proportion of his time within 200 m (mean daily distance = 135 m).
The males with a greater median distance from the nest they sired have limited data (only one day).
We can also summarise the daily min, mean, median and max distance of successful males from nests they sired.
nest_male_date %>%
filter(success==1) %>%
ungroup() %>%
summarise(min_dist=min(dist_m),
max_dist=max(dist_m),
mean_dist=mean(dist_m),
median_dist=median(dist_m),
sd_dist=sd(dist_m)) %>%
flextable()min_dist | max_dist | mean_dist | median_dist | sd_dist | geometry |
|---|---|---|---|---|---|
8.7 [m] | 1,350 [m] | 210 [m] | 116 [m] | 270 | [[XY]] |
In a preliminary analysis (not included here), I explored the three cases where males sired nests far from their territory. P1103 was sired by a male who had a territory further away, but did visit the nest site.
P701 was sired by two males who had neighbouring territories on the other side of the site. The mother was seen on their territories when she was fertile, but apparently decided to build her nest somewhere else.
Determine likely paternity of predated nests
Males with a median daily distance <=200 metres from a nest site still seem to have a higher probability of siring the eggs. We use this information to get a list of “likely” fathers of the predated nests.
potential_fathers <- nest_male_date %>%
filter(nest %in% predated_nests$nest) %>%
group_by(nest, colour_combo) %>%
summarise(first_dist = as.numeric(first(dist_m)),
median_dist = as.numeric(median(dist_m)),
mean_dist = as.numeric(mean(dist_m))) %>%
filter(median_dist<=200) %>%
ungroup() %>%
group_by(colour_combo) %>%
summarise(potential_n_extra_nests=n())
# to get rid of geometry from object
potential_fathers <- data.frame(colour_combo = potential_fathers$colour_combo,
potential_n_extra_nests = potential_fathers$potential_n_extra_nests)
potential_fathers %>% flextable()colour_combo | potential_n_extra_nests |
|---|---|
DB,R,O | 2 |
DB,R,R | 1 |
DG,PI,Y | 1 |
DG,R,DG | 1 |
O,DB,O | 1 |
O,PI,Y | 1 |
O,Y,DG | 4 |
O,Y,O | 2 |
PI,DG,DG | 2 |
PI,O,O | 2 |
PI,PI,PI | 4 |
PI,R,DB | 2 |
PI,Y,R | 3 |
R,DG,R | 1 |
R,O,R | 1 |
R,Y,Y | 2 |
Y,Y,Y | 1 |
write.csv(potential_fathers, file=here("outputs","processed-data","potential_fathers.csv"), row.names=FALSE)potential_fathers <- read.csv(here("outputs","processed-data","potential_fathers.csv"))
father_summary <- full_join(fathers,
potential_fathers,
by=join_by('LL.father'=='colour_combo')) %>%
left_join(.,
tagged_males[,c('colour_combo','bird_ID')],
by=join_by('LL.father'=='colour_combo')) %>%
ungroup() %>%
select(!ID.father) %>%
rename(ID.father=bird_ID) %>%
mutate(success=case_when(n_eggs>0 ~ "yes",
TRUE ~ "no"),
potential_n_extra_nests=replace_na(potential_n_extra_nests,0),
likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
TRUE ~ "less likely"),
likely_undetected_success=case_when(success=="no" & likely_extra_nests=="yes" ~ "yes",
success=="yes" ~ "no",
TRUE ~ "no")) %>%
filter(!is.na(LL.father))Combine success and activity data
activity <- left_join(activity,
father_summary,
by=join_by('bird_ID'=='ID.father',
'colour_combo'=='LL.father')) %>%
left_join(., tagged_males[,c('tag_ID','body_mass')], join_by('tag_ID')) %>%
mutate(n_eggs=replace_na(n_eggs, 0),
n_nests=replace_na(n_nests, 0),
potential_n_extra_nests=replace_na(potential_n_extra_nests,0),
likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
TRUE ~ "less likely"),
success=case_when(n_eggs > 0 ~ "yes",
TRUE ~ "no"),
likely_undetected_success=case_when(success=="no" & likely_extra_nests=="yes" ~ "yes",
TRUE ~ "no"))
daily_activity <- left_join(daily_activity,
father_summary,
by=join_by('bird_ID'=='ID.father',
'colour_combo'=='LL.father')) %>%
left_join(., tagged_males[,c('tag_ID','body_mass')], join_by('tag_ID')) %>%
mutate(n_eggs=replace_na(n_eggs,0),
n_nests=replace_na(n_nests, 0),
potential_n_extra_nests=replace_na(potential_n_extra_nests,0),
likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
TRUE ~ "less likely"),
success=case_when(n_eggs > 0 ~ "yes",
TRUE ~ "no"),
likely_undetected_success=case_when(success=="no" & likely_extra_nests=="yes" ~ "yes",
TRUE ~ "no"))activity <- activity %>%
mutate(breeding_stage=case_when(window_end_localtime-10*60<=start_FertilePeriod ~ "pre-fertile",
window_end_localtime-10*60>start_FertilePeriod &
window_end_localtime<end_FertilePeriod ~ "fertile",
window_end_localtime>=end_FertilePeriod ~ "post-fertile"))
daily_activity <- daily_activity %>%
mutate(breeding_stage=case_when(local_date<=start_FertilePeriod ~ "pre-fertile",
local_date>start_FertilePeriod &
local_date<end_FertilePeriod ~ "fertile",
local_date>=end_FertilePeriod ~ "post-fertile"))
daily_activity <- left_join(daily_activity,
ff_by_date,
by=join_by('local_date'=='date'),
relationship = "many-to-one") %>%
mutate(n_FertileFemales=replace_na(n_FertileFemales,0))Summarise individual data and exclude incomplete data
daily_activity <- daily_activity %>%
filter(n_hours>=min_acc_hrs_per_day)ind_total_active <- daily_activity %>%
group_by(bird_ID, colour_combo, breeding_stage,
body_mass, n_eggs, n_nests, success,
likely_undetected_success, potential_n_extra_nests, likely_extra_nests) %>%
summarise(n_windows=sum(n_windows),
n_days=n(),
total_perc_active=sum(total_active)/sum(n_windows)*100)
# Note: total perc active will be same whether calculated from activity or daily_activity
# However, mean ODBA should be calculated from activity
ind_odba <- activity %>%
filter(!is.na(activity_threshold)) %>%
group_by(bird_ID, colour_combo, breeding_stage,
body_mass, n_eggs, n_nests, success,
likely_undetected_success, potential_n_extra_nests, likely_extra_nests) %>%
summarise(mean_ODBA=mean(mean_ODBA))
ind_activity <- full_join(ind_total_active, ind_odba,
by=c('bird_ID','colour_combo','breeding_stage',
'body_mass','n_eggs','n_nests','success',
'likely_undetected_success','potential_n_extra_nests',
'likely_extra_nests'))
# Filter out missing/insufficient data
ind_activity <- ind_activity %>% filter(!is.na(total_perc_active))
rm(ind_total_active, ind_odba)Note: There are no males with <2 days of data during the fertile period and >=2 days of data during post-fertile period. None of the males with <2 days of data during the fertile period were successful.
males_with_sufficient_data = ind_activity %>%
filter(breeding_stage=="fertile" & n_days>=min_acc_days) %>%
select(colour_combo, bird_ID, n_days, n_windows, success)
males_with_sufficient_data %>%
arrange(n_days, n_windows, success) %>%
select(colour_combo, bird_ID, n_days, n_windows, n_nests) %>%
flextable()breeding_stage | body_mass | n_eggs | success | likely_undetected_success | potential_n_extra_nests | colour_combo | bird_ID | n_days | n_windows | n_nests |
|---|---|---|---|---|---|---|---|---|---|---|
fertile | 99.5 | 0 | no | no | 0 | DB,PI,DB | 141,269,115 | 2 | 268 | 0 |
fertile | 97.4 | 0 | no | no | 0 | Y,PI,DB | 141,269,180 | 2 | 275 | 0 |
fertile | 95.4 | 0 | no | no | 0 | PI,O,PI | 141,269,158 | 2 | 280 | 0 |
fertile | 105.6 | 0 | no | no | 0 | DG,PI,DB | 141,269,037 | 2 | 288 | 0 |
fertile | 105.7 | 0 | no | no | 0 | R,R,Y | 141,269,132 | 2 | 288 | 0 |
fertile | 106.7 | 0 | no | no | 0 | R,DB,O | 141,269,039 | 3 | 396 | 0 |
fertile | 90.4 | 0 | no | no | 0 | R,Y,DB | 141,269,159 | 3 | 417 | 0 |
fertile | 106.4 | 0 | no | no | 0 | DB,DB,O | 141,269,031 | 3 | 431 | 0 |
fertile | 109.8 | 0 | no | no | 0 | O,R,Y | 141,269,128 | 3 | 432 | 0 |
fertile | 101.4 | 0 | no | yes | 2 | R,Y,Y | 141,269,173 | 3 | 432 | 0 |
fertile | 96.6 | 4 | yes | no | 0 | DG,DB,DG | 141,269,127 | 3 | 432 | 1 |
fertile | 93.7 | 0 | no | no | 0 | PI,Y,PI | 141,269,145 | 4 | 553 | 0 |
fertile | 93.5 | 0 | no | no | 0 | DG,Y,Y | 141,269,038 | 4 | 569 | 0 |
fertile | 96.7 | 0 | no | no | 0 | DB,DB,Y | 141,269,179 | 4 | 570 | 0 |
fertile | 101.5 | 0 | no | no | 0 | O,DG,R | 141,269,004 | 4 | 576 | 0 |
fertile | 99.3 | 0 | no | no | 0 | DG,DB,PI | 141,269,163 | 4 | 576 | 0 |
fertile | 89.0 | 0 | no | no | 0 | DG,PI,O | 141,269,050 | 5 | 712 | 0 |
fertile | 108.2 | 0 | no | no | 0 | R,DG,Y | 141,269,048 | 5 | 718 | 0 |
fertile | 108.3 | 0 | no | no | 0 | O,R,PI | 141,269,047 | 5 | 720 | 0 |
fertile | 105.5 | 0 | no | no | 0 | Y,R,PI | 141,269,151 | 5 | 720 | 0 |
fertile | 103.5 | 0 | no | yes | 4 | O,Y,DG | 141,269,155 | 5 | 720 | 0 |
fertile | 99.2 | 0 | no | no | 0 | PI,DB,DB | 141,269,183 | 5 | 720 | 0 |
fertile | 95.3 | 4 | yes | no | 0 | DB,DG,Y | 141,269,046 | 5 | 720 | 1 |
fertile | 95.2 | 0 | no | no | 0 | R,O,O | 141,269,045 | 6 | 843 | 0 |
fertile | 97.9 | 2 | yes | no | 2 | O,Y,O | 141,269,182 | 6 | 845 | 1 |
fertile | 92.5 | 5 | yes | no | 0 | Y,DG,PI | 141,269,044 | 6 | 846 | 2 |
fertile | 93.2 | 0 | no | no | 0 | R,Y,R | 141,269,181 | 6 | 848 | 0 |
fertile | 101.1 | 4 | yes | no | 0 | DG,O,O | 141,269,043 | 6 | 857 | 1 |
fertile | 103.5 | 0 | no | no | 0 | PI,PI,DB | 141,269,178 | 6 | 861 | 0 |
fertile | 105.7 | 0 | no | no | 0 | R,Y,O | 141,269,042 | 6 | 863 | 0 |
fertile | 97.4 | 0 | no | no | 0 | Y,Y,O | 141,269,041 | 6 | 864 | 0 |
fertile | 102.5 | 0 | no | no | 0 | Y,DG,DG | 141,269,167 | 6 | 864 | 0 |
fertile | 97.6 | 0 | no | no | 0 | DB,PI,PI | 141,269,171 | 6 | 864 | 0 |
fertile | 111.4 | 0 | no | yes | 1 | DG,PI,Y | 141,269,172 | 6 | 864 | 0 |
fertile | 87.7 | 0 | no | yes | 4 | PI,PI,PI | 141,269,174 | 6 | 864 | 0 |
fertile | 105.4 | 0 | no | no | 0 | DB,R,Y | 141,269,177 | 6 | 864 | 0 |
fertile | 98.5 | 3 | yes | no | 0 | DG,R,R | 141,269,175 | 6 | 864 | 1 |
fertile | 99.9 | 0 | no | no | 0 | Y,Y,PI | 141,269,040 | 7 | 987 | 0 |
fertile | 108.7 | 0 | no | no | 0 | R,O,DB | 141,269,136 | 7 | 987 | 0 |
fertile | 96.7 | 0 | no | no | 0 | R,DB,R | 141,269,170 | 7 | 989 | 0 |
fertile | 97.5 | 4 | yes | no | 1 | DB,R,R | 141,269,169 | 7 | 1,000 | 1 |
fertile | 97.9 | 0 | no | no | 0 | DB,DG,DG | 141,269,165 | 7 | 1,007 | 0 |
fertile | 93.8 | 0 | no | no | 0 | DB,DB,DB | 141,269,006 | 7 | 1,008 | 0 |
fertile | 95.7 | 0 | no | no | 0 | DG,R,DB | 141,269,034 | 7 | 1,008 | 0 |
fertile | 102.5 | 0 | no | no | 0 | R,PI,PI | 141,269,036 | 7 | 1,008 | 0 |
fertile | 102.4 | 0 | no | no | 0 | DG,DG,DB | 141,269,133 | 7 | 1,008 | 0 |
fertile | 103.5 | 0 | no | no | 0 | O,DB,R | 141,269,160 | 7 | 1,008 | 0 |
fertile | 115.9 | 0 | no | no | 0 | DG,DB,Y | 141,269,162 | 7 | 1,008 | 0 |
fertile | 94.3 | 0 | no | no | 0 | O,PI,O | 141,269,164 | 7 | 1,008 | 0 |
fertile | 100.4 | 0 | no | yes | 2 | PI,O,O | 141,269,166 | 7 | 1,008 | 0 |
fertile | 107.3 | 0 | no | no | 0 | DG,DG,R | 141,269,033 | 8 | 1,144 | 0 |
fertile | 104.7 | 0 | no | yes | 2 | PI,R,DB | 141,269,157 | 8 | 1,146 | 0 |
fertile | 93.2 | 0 | no | yes | 1 | O,PI,Y | 141,269,156 | 8 | 1,150 | 0 |
fertile | 103.9 | 0 | no | no | 0 | DG,O,DG | 141,269,030 | 8 | 1,152 | 0 |
fertile | 105.5 | 0 | no | yes | 1 | DG,R,DG | 141,269,149 | 8 | 1,152 | 0 |
fertile | 99.5 | 0 | no | no | 0 | DG,R,O | 141,269,153 | 8 | 1,152 | 0 |
fertile | 100.4 | 0 | no | no | 0 | O,R,O | 141,269,154 | 8 | 1,152 | 0 |
fertile | 108.1 | 3 | yes | no | 2 | DB,R,O | 141,269,148 | 8 | 1,152 | 1 |
fertile | 100.5 | 1 | yes | no | 1 | Y,Y,Y | 141,269,152 | 8 | 1,152 | 1 |
fertile | 105.2 | 5 | yes | no | 0 | O,DG,DG | 141,269,147 | 9 | 1,276 | 2 |
fertile | 101.8 | 2 | yes | no | 0 | R,DG,PI | 141,269,146 | 9 | 1,281 | 1 |
fertile | 102.4 | 3 | yes | no | 0 | DB,Y,R | 141,269,144 | 9 | 1,289 | 1 |
fertile | 102.2 | 0 | no | no | 0 | Y,O,R | 141,269,143 | 9 | 1,292 | 0 |
fertile | 106.9 | 0 | no | yes | 1 | R,DG,R | 141,269,137 | 9 | 1,295 | 0 |
fertile | 97.5 | 0 | no | yes | 1 | R,O,R | 141,269,138 | 9 | 1,296 | 0 |
fertile | 89.2 | 0 | no | no | 0 | PI,O,R | 141,269,140 | 9 | 1,296 | 0 |
fertile | 88.4 | 0 | no | no | 0 | R,R,DG | 141,269,141 | 9 | 1,296 | 0 |
fertile | 94.7 | 4 | yes | no | 1 | O,DB,O | 141,269,026 | 9 | 1,296 | 1 |
fertile | 97.2 | 0 | no | yes | 3 | PI,Y,R | 141,269,011 | 10 | 1,422 | 0 |
fertile | 99.6 | 0 | no | no | 0 | DB,Y,PI | 141,269,134 | 10 | 1,436 | 0 |
fertile | 109.6 | 0 | no | no | 0 | PI,PI,O | 141,269,024 | 10 | 1,437 | 0 |
fertile | 97.5 | 0 | no | no | 0 | Y,DB,Y | 141,269,022 | 10 | 1,440 | 0 |
fertile | 104.3 | 0 | no | no | 0 | Y,Y,DG | 141,269,129 | 10 | 1,440 | 0 |
fertile | 103.2 | 0 | no | no | 0 | Y,O,Y | 141,269,126 | 11 | 1,566 | 0 |
fertile | 99.2 | 0 | no | yes | 2 | PI,DG,DG | 141,269,122 | 11 | 1,584 | 0 |
fertile | 95.8 | 1 | yes | no | 0 | Y,Y,R | 141,269,125 | 11 | 1,584 | 1 |
fertile | 93.5 | 0 | no | no | 0 | PI,Y,Y | 141,269,105 | 12 | 1,720 | 0 |
fertile | 100.1 | 0 | no | no | 0 | O,DB,DG | 141,269,010 | 12 | 1,728 | 0 |
fertile | 107.5 | 0 | no | no | 0 | DB,O,DG | 141,269,118 | 12 | 1,728 | 0 |
fertile | 101.4 | 6 | yes | no | 0 | PI,PI,R | 141,269,009 | 12 | 1,728 | 2 |
fertile | 106.0 | 3 | yes | no | 0 | PI,R,PI | 141,269,114 | 13 | 1,871 | 1 |
There are 81 males with at least 2 days of complete (20.4 hours) data during the fertile period, as shown in the table above.
From here on, the activity datasets are restricted only to these males.
daily_activity <- daily_activity %>%
filter(bird_ID %in% males_with_sufficient_data$bird_ID &
n_hours>=min_acc_hrs_per_day)
ind_activity <- ind_activity %>%
filter(bird_ID %in% males_with_sufficient_data$bird_ID)ind_activity <- left_join(ind_activity,
male_tenure,
join_by('bird_ID', 'colour_combo'))Q: Do males spend more time active when there are more mating opportunities?
The percentage of time that males spent active in a day ranged from less than 20% to over 99% (mean ± standard deviation: 55.634 +/- 21.494.
daily_activity %>%
ungroup() %>%
summarise(min_perc_active=min(perc_active),
max_perc_active=max(perc_active),
mean_perc_active=mean(perc_active),
sd_perc_active=sd(perc_active),
n_birds=n_distinct(bird_ID),
n=n(),
se_perc_active=sd_perc_active/sqrt(n)) %>%
flextable()min_perc_active | max_perc_active | mean_perc_active | sd_perc_active | n_birds | n | se_perc_active |
|---|---|---|---|---|---|---|
4.86 | 100 | 55.6 | 21.5 | 81 | 802 | 0.759 |
within_ind_var <- daily_activity %>%
ungroup() %>%
select(colour_combo, perc_active, mean_ODBA) %>%
group_by(colour_combo) %>%
summarise(min_perc_active=min(perc_active),
max_perc_active=max(perc_active),
diff_perc_active=max_perc_active-min_perc_active,
diff_hours_active=diff_perc_active*24/100,
diff_multiply_active=max_perc_active/min_perc_active,
mean_perc_active=mean(perc_active),
sd_perc_active=sd(perc_active),
min_ODBA=min(mean_ODBA),
max_ODBA=max(mean_ODBA),
diff_ODBA=max_ODBA-min_ODBA,
diff_multiply_ODBA=max_ODBA/min_ODBA,
mean_mean_ODBA=mean(mean_ODBA),
sd_ODBA=sd(mean_ODBA),
n_days=n()) %>%
arrange(diff_multiply_ODBA)
within_ind_var %>% flextable()colour_combo | min_perc_active | max_perc_active | diff_perc_active | diff_hours_active | diff_multiply_active | mean_perc_active | sd_perc_active | min_ODBA | max_ODBA | diff_ODBA | diff_multiply_ODBA | mean_mean_ODBA | sd_ODBA | n_days |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Y,PI,DB | 55.73 | 63.2 | 7.47 | 1.793 | 1.13 | 59.5 | 5.28 | 0.2369 | 0.250 | 0.0129 | 1.05 | 0.243 | 0.00913 | 2 |
DB,PI,DB | 58.06 | 60.4 | 2.35 | 0.565 | 1.04 | 59.2 | 1.66 | 0.3175 | 0.347 | 0.0292 | 1.09 | 0.332 | 0.02061 | 2 |
O,DG,R | 63.19 | 75.0 | 11.81 | 2.833 | 1.19 | 69.6 | 5.93 | 0.3069 | 0.345 | 0.0377 | 1.12 | 0.333 | 0.01768 | 4 |
DB,DB,DB | 49.31 | 63.9 | 14.58 | 3.500 | 1.30 | 58.1 | 6.12 | 0.2453 | 0.308 | 0.0631 | 1.26 | 0.283 | 0.02566 | 7 |
PI,Y,PI | 48.09 | 77.1 | 28.99 | 6.958 | 1.60 | 60.5 | 12.14 | 0.2521 | 0.323 | 0.0708 | 1.28 | 0.280 | 0.03024 | 4 |
O,PI,O | 39.58 | 54.9 | 15.28 | 3.667 | 1.39 | 48.5 | 5.68 | 0.1944 | 0.251 | 0.0564 | 1.29 | 0.222 | 0.02284 | 7 |
R,Y,DB | 31.25 | 63.9 | 32.64 | 7.833 | 2.04 | 48.3 | 16.36 | 0.2769 | 0.359 | 0.0820 | 1.30 | 0.305 | 0.04651 | 3 |
DG,PI,DB | 38.19 | 54.2 | 15.97 | 3.833 | 1.42 | 46.2 | 11.29 | 0.1918 | 0.249 | 0.0576 | 1.30 | 0.221 | 0.04075 | 2 |
Y,R,PI | 41.67 | 62.5 | 20.83 | 5.000 | 1.50 | 53.2 | 7.63 | 0.1424 | 0.186 | 0.0436 | 1.31 | 0.165 | 0.01587 | 5 |
Y,O,Y | 27.08 | 52.8 | 25.69 | 6.167 | 1.95 | 41.1 | 6.96 | 0.1656 | 0.217 | 0.0517 | 1.31 | 0.189 | 0.01862 | 13 |
Y,Y,O | 40.28 | 61.1 | 20.83 | 5.000 | 1.52 | 51.5 | 6.80 | 0.1654 | 0.223 | 0.0572 | 1.35 | 0.198 | 0.01898 | 10 |
DG,DB,DG | 36.11 | 59.0 | 22.92 | 5.500 | 1.63 | 49.3 | 11.85 | 0.1972 | 0.268 | 0.0708 | 1.36 | 0.240 | 0.03756 | 3 |
DB,DB,O | 46.53 | 67.1 | 20.61 | 4.945 | 1.44 | 57.8 | 10.44 | 0.2191 | 0.303 | 0.0840 | 1.38 | 0.256 | 0.04302 | 3 |
DB,Y,PI | 29.86 | 47.9 | 18.06 | 4.333 | 1.60 | 39.4 | 6.42 | 0.1486 | 0.212 | 0.0633 | 1.43 | 0.178 | 0.02132 | 12 |
DG,R,DB | 40.97 | 66.7 | 25.69 | 6.167 | 1.63 | 56.2 | 6.98 | 0.1841 | 0.268 | 0.0839 | 1.46 | 0.231 | 0.02551 | 9 |
R,DG,Y | 36.62 | 59.0 | 22.41 | 5.378 | 1.61 | 49.6 | 8.20 | 0.1663 | 0.243 | 0.0763 | 1.46 | 0.210 | 0.02835 | 10 |
O,R,Y | 52.78 | 64.6 | 11.81 | 2.833 | 1.22 | 60.4 | 6.62 | 0.2075 | 0.303 | 0.0954 | 1.46 | 0.257 | 0.04783 | 3 |
R,R,Y | 27.78 | 51.4 | 23.61 | 5.667 | 1.85 | 39.6 | 16.70 | 0.1542 | 0.232 | 0.0782 | 1.51 | 0.193 | 0.05529 | 2 |
DB,DB,Y | 58.33 | 91.3 | 32.97 | 7.913 | 1.57 | 73.7 | 14.02 | 0.3353 | 0.512 | 0.1766 | 1.53 | 0.417 | 0.07697 | 4 |
DG,DB,Y | 38.19 | 65.3 | 27.08 | 6.500 | 1.71 | 49.7 | 10.69 | 0.1642 | 0.253 | 0.0892 | 1.54 | 0.200 | 0.03196 | 9 |
Y,DG,DG | 31.25 | 56.9 | 25.69 | 6.167 | 1.82 | 47.3 | 8.72 | 0.1509 | 0.234 | 0.0835 | 1.55 | 0.207 | 0.02755 | 7 |
DG,DB,PI | 39.58 | 70.1 | 30.56 | 7.333 | 1.77 | 56.4 | 13.71 | 0.1915 | 0.300 | 0.1086 | 1.57 | 0.250 | 0.05316 | 4 |
DG,R,O | 26.39 | 50.0 | 23.61 | 5.667 | 1.89 | 39.1 | 8.33 | 0.1194 | 0.188 | 0.0681 | 1.57 | 0.146 | 0.02204 | 11 |
DB,R,Y | 23.61 | 54.2 | 30.56 | 7.333 | 2.29 | 37.3 | 9.10 | 0.1300 | 0.206 | 0.0758 | 1.58 | 0.162 | 0.02303 | 9 |
PI,O,PI | 4.86 | 26.5 | 21.61 | 5.186 | 5.45 | 15.7 | 15.28 | 0.1356 | 0.217 | 0.0812 | 1.60 | 0.176 | 0.05741 | 2 |
DG,PI,O | 31.94 | 53.4 | 21.44 | 5.145 | 1.67 | 43.9 | 8.02 | 0.1389 | 0.224 | 0.0851 | 1.61 | 0.177 | 0.02543 | 9 |
PI,PI,DB | 14.18 | 39.6 | 25.40 | 6.096 | 2.79 | 34.5 | 10.05 | 0.1099 | 0.179 | 0.0688 | 1.63 | 0.160 | 0.02548 | 6 |
R,Y,Y | 39.58 | 73.6 | 34.03 | 8.167 | 1.86 | 52.1 | 11.31 | 0.1779 | 0.294 | 0.1157 | 1.65 | 0.225 | 0.04828 | 7 |
R,DB,O | 17.83 | 49.6 | 31.76 | 7.623 | 2.78 | 31.0 | 16.54 | 0.1331 | 0.221 | 0.0874 | 1.66 | 0.174 | 0.04393 | 3 |
R,O,DB | 19.51 | 48.6 | 29.10 | 6.984 | 2.49 | 39.4 | 10.29 | 0.1099 | 0.183 | 0.0733 | 1.67 | 0.162 | 0.02591 | 7 |
DG,Y,Y | 42.66 | 75.0 | 32.34 | 7.762 | 1.76 | 57.8 | 13.43 | 0.2381 | 0.400 | 0.1614 | 1.68 | 0.314 | 0.06655 | 4 |
DB,PI,PI | 40.28 | 75.0 | 34.72 | 8.333 | 1.86 | 55.6 | 11.11 | 0.1386 | 0.241 | 0.1020 | 1.74 | 0.181 | 0.02787 | 14 |
O,PI,Y | 22.92 | 58.3 | 35.42 | 8.500 | 2.55 | 39.9 | 10.29 | 0.1375 | 0.239 | 0.1017 | 1.74 | 0.181 | 0.02529 | 24 |
PI,DB,DB | 30.37 | 75.7 | 45.32 | 10.878 | 2.49 | 50.6 | 15.24 | 0.1575 | 0.275 | 0.1173 | 1.74 | 0.204 | 0.03921 | 10 |
O,R,PI | 44.44 | 77.8 | 33.33 | 8.000 | 1.75 | 62.4 | 12.27 | 0.1922 | 0.336 | 0.1437 | 1.75 | 0.260 | 0.04871 | 8 |
DG,R,R | 31.94 | 72.2 | 40.28 | 9.667 | 2.26 | 52.2 | 13.64 | 0.2296 | 0.402 | 0.1727 | 1.75 | 0.313 | 0.06048 | 8 |
PI,O,R | 36.11 | 72.9 | 36.81 | 8.833 | 2.02 | 50.6 | 12.38 | 0.1767 | 0.328 | 0.1514 | 1.86 | 0.228 | 0.05180 | 11 |
PI,Y,R | 11.11 | 59.0 | 47.92 | 11.500 | 5.31 | 37.7 | 11.51 | 0.1314 | 0.245 | 0.1140 | 1.87 | 0.193 | 0.02576 | 20 |
Y,DG,PI | 29.17 | 68.8 | 39.58 | 9.500 | 2.36 | 48.5 | 11.69 | 0.1375 | 0.258 | 0.1209 | 1.88 | 0.193 | 0.03527 | 16 |
DG,DG,DB | 52.78 | 98.6 | 45.83 | 11.000 | 1.87 | 78.0 | 19.35 | 0.3233 | 0.682 | 0.3584 | 2.11 | 0.498 | 0.15002 | 7 |
DG,R,DG | 19.44 | 68.8 | 49.31 | 11.833 | 3.54 | 51.0 | 13.63 | 0.1174 | 0.254 | 0.1365 | 2.16 | 0.202 | 0.03907 | 17 |
O,Y,DG | 66.67 | 99.3 | 32.64 | 7.833 | 1.49 | 81.5 | 12.41 | 0.3059 | 0.663 | 0.3574 | 2.17 | 0.432 | 0.14224 | 5 |
Y,Y,Y | 58.33 | 99.3 | 40.97 | 9.833 | 1.70 | 79.4 | 15.48 | 0.2734 | 0.602 | 0.3288 | 2.20 | 0.409 | 0.12318 | 10 |
O,DB,O | 37.50 | 91.7 | 54.17 | 13.000 | 2.44 | 63.4 | 18.19 | 0.2571 | 0.574 | 0.3171 | 2.23 | 0.381 | 0.09446 | 11 |
R,R,DG | 50.69 | 100.0 | 49.31 | 11.833 | 1.97 | 73.7 | 15.03 | 0.2520 | 0.583 | 0.3306 | 2.31 | 0.371 | 0.09908 | 11 |
PI,O,O | 36.81 | 97.9 | 61.11 | 14.667 | 2.66 | 67.4 | 21.93 | 0.2311 | 0.566 | 0.3351 | 2.45 | 0.353 | 0.12218 | 7 |
DB,Y,R | 31.94 | 91.7 | 59.72 | 14.333 | 2.87 | 57.1 | 18.44 | 0.1691 | 0.414 | 0.2454 | 2.45 | 0.258 | 0.08187 | 12 |
PI,R,DB | 46.53 | 100.0 | 53.47 | 12.833 | 2.15 | 73.7 | 17.86 | 0.2597 | 0.642 | 0.3826 | 2.47 | 0.409 | 0.13409 | 9 |
R,DB,R | 34.03 | 84.0 | 50.00 | 12.000 | 2.47 | 57.3 | 12.46 | 0.1640 | 0.413 | 0.2491 | 2.52 | 0.252 | 0.06536 | 15 |
DG,O,O | 35.42 | 83.9 | 48.52 | 11.646 | 2.37 | 59.4 | 14.54 | 0.1905 | 0.500 | 0.3092 | 2.62 | 0.294 | 0.08758 | 10 |
DB,DG,DG | 20.14 | 90.2 | 70.07 | 16.817 | 4.48 | 49.9 | 20.88 | 0.1682 | 0.441 | 0.2733 | 2.62 | 0.282 | 0.08086 | 10 |
R,Y,O | 33.33 | 93.0 | 59.67 | 14.322 | 2.79 | 55.8 | 23.18 | 0.1753 | 0.462 | 0.2869 | 2.64 | 0.270 | 0.11065 | 6 |
R,PI,PI | 32.64 | 96.5 | 63.89 | 15.333 | 2.96 | 63.0 | 22.96 | 0.2363 | 0.634 | 0.3977 | 2.68 | 0.376 | 0.13598 | 8 |
PI,R,PI | 18.75 | 89.6 | 70.83 | 17.000 | 4.78 | 66.0 | 18.05 | 0.2153 | 0.580 | 0.3645 | 2.69 | 0.400 | 0.08569 | 16 |
PI,PI,O | 36.11 | 98.6 | 62.50 | 15.000 | 2.73 | 70.0 | 18.53 | 0.2158 | 0.581 | 0.3654 | 2.69 | 0.369 | 0.10271 | 14 |
O,DB,R | 9.03 | 71.5 | 62.50 | 15.000 | 7.92 | 47.8 | 17.05 | 0.1008 | 0.271 | 0.1707 | 2.69 | 0.197 | 0.05015 | 11 |
Y,Y,R | 22.92 | 94.4 | 71.53 | 17.167 | 4.12 | 61.9 | 25.87 | 0.1849 | 0.498 | 0.3131 | 2.69 | 0.320 | 0.12107 | 13 |
R,DG,PI | 39.58 | 98.6 | 59.03 | 14.167 | 2.49 | 70.1 | 17.17 | 0.2219 | 0.599 | 0.3774 | 2.70 | 0.368 | 0.11371 | 10 |
PI,PI,R | 5.56 | 100.0 | 94.44 | 22.667 | 18.00 | 50.4 | 31.79 | 0.2097 | 0.604 | 0.3944 | 2.88 | 0.391 | 0.13949 | 15 |
DG,O,DG | 30.56 | 97.2 | 66.67 | 16.000 | 3.18 | 63.8 | 23.44 | 0.1700 | 0.495 | 0.3248 | 2.91 | 0.301 | 0.10971 | 11 |
Y,O,R | 15.97 | 92.4 | 76.39 | 18.333 | 5.78 | 64.4 | 22.35 | 0.1753 | 0.513 | 0.3372 | 2.92 | 0.363 | 0.09593 | 14 |
PI,Y,Y | 31.62 | 98.6 | 66.99 | 16.078 | 3.12 | 76.1 | 19.28 | 0.2359 | 0.699 | 0.4632 | 2.96 | 0.489 | 0.13686 | 12 |
DB,O,DG | 7.64 | 95.8 | 88.19 | 21.167 | 12.55 | 44.7 | 25.67 | 0.1652 | 0.506 | 0.3409 | 3.06 | 0.293 | 0.10755 | 22 |
DB,DG,Y | 14.58 | 79.2 | 64.58 | 15.500 | 5.43 | 35.8 | 20.08 | 0.1253 | 0.390 | 0.2649 | 3.11 | 0.191 | 0.08119 | 10 |
O,DB,DG | 17.36 | 98.6 | 81.25 | 19.500 | 5.68 | 62.6 | 26.71 | 0.1691 | 0.536 | 0.3670 | 3.17 | 0.348 | 0.12297 | 22 |
Y,DB,Y | 18.75 | 92.4 | 73.61 | 17.667 | 4.93 | 53.5 | 25.28 | 0.1608 | 0.516 | 0.3553 | 3.21 | 0.284 | 0.12621 | 18 |
R,O,O | 22.22 | 88.6 | 66.40 | 15.935 | 3.99 | 43.0 | 22.17 | 0.1052 | 0.341 | 0.2355 | 3.24 | 0.178 | 0.07560 | 9 |
O,Y,O | 25.69 | 88.0 | 62.31 | 14.953 | 3.42 | 61.7 | 24.14 | 0.1528 | 0.498 | 0.3455 | 3.26 | 0.334 | 0.13877 | 8 |
O,DG,DG | 25.00 | 98.6 | 73.60 | 17.664 | 3.94 | 73.9 | 21.22 | 0.1613 | 0.528 | 0.3664 | 3.27 | 0.364 | 0.10335 | 15 |
DB,R,R | 38.89 | 97.8 | 58.91 | 14.137 | 2.51 | 67.8 | 20.09 | 0.1943 | 0.642 | 0.4478 | 3.30 | 0.341 | 0.12866 | 11 |
PI,DG,DG | 31.94 | 98.6 | 66.67 | 16.000 | 3.09 | 64.5 | 24.22 | 0.1797 | 0.599 | 0.4192 | 3.33 | 0.350 | 0.15376 | 16 |
DG,PI,Y | 36.11 | 98.6 | 62.50 | 15.000 | 2.73 | 61.5 | 21.20 | 0.1895 | 0.637 | 0.4472 | 3.36 | 0.323 | 0.14337 | 10 |
PI,PI,PI | 24.31 | 96.5 | 72.22 | 17.333 | 3.97 | 67.2 | 19.11 | 0.1601 | 0.552 | 0.3918 | 3.45 | 0.339 | 0.09913 | 11 |
O,R,O | 34.72 | 98.6 | 63.89 | 15.333 | 2.84 | 65.3 | 23.69 | 0.1588 | 0.548 | 0.3891 | 3.45 | 0.309 | 0.14162 | 10 |
DG,DG,R | 24.31 | 97.9 | 73.61 | 17.667 | 4.03 | 58.4 | 26.36 | 0.1588 | 0.557 | 0.3979 | 3.51 | 0.305 | 0.13729 | 12 |
R,O,R | 37.50 | 97.9 | 60.42 | 14.500 | 2.61 | 66.7 | 21.59 | 0.1706 | 0.601 | 0.4301 | 3.52 | 0.334 | 0.14649 | 12 |
Y,Y,PI | 25.00 | 96.0 | 70.97 | 17.032 | 3.84 | 50.8 | 25.68 | 0.1199 | 0.428 | 0.3077 | 3.57 | 0.225 | 0.12562 | 8 |
Y,Y,DG | 15.97 | 97.9 | 81.94 | 19.667 | 6.13 | 44.3 | 27.91 | 0.1535 | 0.655 | 0.5020 | 4.27 | 0.277 | 0.15214 | 20 |
DB,R,O | 11.81 | 99.3 | 87.50 | 21.000 | 8.41 | 61.6 | 29.07 | 0.1249 | 0.578 | 0.4536 | 4.63 | 0.316 | 0.15355 | 14 |
R,Y,R | 14.58 | 90.3 | 75.69 | 18.167 | 6.19 | 49.2 | 28.02 | 0.0967 | 0.505 | 0.4078 | 5.22 | 0.245 | 0.14058 | 14 |
R,DG,R | 19.44 | 100.0 | 80.56 | 19.333 | 5.14 | 60.0 | 31.50 | 0.1416 | 0.807 | 0.6657 | 5.70 | 0.316 | 0.20293 | 12 |
daily_activity %>%
ungroup() %>%
group_by(breeding_stage) %>%
summarise(mean_perc_active=mean(perc_active),
sd_perc_active=sd(perc_active),
n_samples=n(),
n_birds=n_distinct(bird_ID),
se_perc_active=sd_perc_active/sqrt(n_samples)) %>%
flextable()breeding_stage | mean_perc_active | sd_perc_active | n_samples | n_birds | se_perc_active |
|---|---|---|---|---|---|
fertile | 62.1 | 20.7 | 555 | 81 | 0.877 |
post-fertile | 41.0 | 15.3 | 247 | 56 | 0.972 |
daily_activity %>%
ungroup() %>%
group_by(breeding_stage) %>%
summarise(min_mean_ODBA=min(mean_ODBA),
max_mean_ODBA=max(mean_ODBA),
mean_mean_ODBA=mean(mean_ODBA),
sd_mean_ODBA=sd(mean_ODBA),
n_samples=n(),
n_birds=n_distinct(bird_ID),
se_mean_ODBA=sd_mean_ODBA/sqrt(n_samples)) %>%
flextable()breeding_stage | min_mean_ODBA | max_mean_ODBA | mean_mean_ODBA | sd_mean_ODBA | n_samples | n_birds | se_mean_ODBA |
|---|---|---|---|---|---|---|---|
fertile | 0.1008 | 0.807 | 0.315 | 0.1320 | 555 | 81 | 0.0056 |
post-fertile | 0.0967 | 0.467 | 0.205 | 0.0518 | 247 | 56 | 0.0033 |
daily_activity %>%
ungroup() %>%
select(colour_combo, perc_active, mean_ODBA, breeding_stage) %>%
group_by(colour_combo, breeding_stage) %>%
summarise(mean_perc_active=mean(perc_active),
mean_ODBA=mean(mean_ODBA)) %>%
pivot_wider(values_from = c('mean_perc_active','mean_ODBA'),
names_from = 'breeding_stage') %>%
filter(!is.na(`mean_ODBA_post-fertile`) & !is.na(mean_perc_active_fertile)) %>%
mutate(diff_perc_active=mean_perc_active_fertile-`mean_perc_active_post-fertile`,
diff_hours_active=24/100*diff_perc_active,
diff_mean_ODBA=mean_ODBA_fertile-`mean_ODBA_post-fertile`) %>%
ungroup() %>%
summarise(mean_diff_perc_active=mean(diff_perc_active),
sd_diff_perc_active=sd(diff_perc_active),
mean_diff_hours_active=mean(diff_hours_active),
sd_diff_hours_active=sd(diff_hours_active),
mean_diff_ODBA=mean(diff_mean_ODBA),
sd_diff_ODBA=sd(diff_mean_ODBA)) %>%
flextable()mean_diff_perc_active | sd_diff_perc_active | mean_diff_hours_active | sd_diff_hours_active | mean_diff_ODBA | sd_diff_ODBA |
|---|---|---|---|---|---|
20.6 | 18.1 | 4.95 | 4.35 | 0.0991 | 0.0832 |
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
select(colour_combo, perc_active, n_FertileFemales) %>%
group_by(colour_combo, n_FertileFemales) %>%
summarise(mean_perc_active=mean(perc_active)) %>%
pivot_wider(values_from = 'mean_perc_active', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
filter(!is.na(ff_40) & !is.na(ff_0)) %>%
mutate(diff=ff_40-ff_0,
diff_hours=24/100*diff) %>%
flextable()colour_combo | ff_40 | ff_0 | diff | diff_hours |
|---|---|---|---|---|
DB,O,DG | 81.2 | 20.0 | 61.23 | 14.694 |
DB,PI,PI | 54.9 | 62.8 | -7.99 | -1.917 |
DB,R,O | 99.0 | 19.4 | 79.51 | 19.083 |
DG,R,DG | 53.8 | 49.9 | 3.90 | 0.936 |
O,DB,DG | 95.5 | 28.8 | 66.67 | 16.000 |
O,DG,DG | 97.9 | 28.8 | 69.09 | 16.582 |
O,PI,Y | 35.9 | 39.1 | -3.20 | -0.768 |
PI,DG,DG | 98.6 | 34.0 | 64.58 | 15.500 |
PI,PI,PI | 96.5 | 24.3 | 72.22 | 17.333 |
PI,Y,R | 38.2 | 25.6 | 12.62 | 3.028 |
R,DB,R | 78.0 | 56.9 | 21.07 | 5.057 |
R,Y,R | 66.4 | 22.1 | 44.27 | 10.624 |
Y,DB,Y | 85.1 | 30.0 | 55.03 | 13.208 |
Y,DG,PI | 57.9 | 46.1 | 11.87 | 2.849 |
Y,O,R | 89.6 | 16.0 | 73.61 | 17.667 |
Y,Y,DG | 69.1 | 22.9 | 46.18 | 11.083 |
The table above shows the percentage of time that males (each with a unique colour_combo) were active when female fertility was at its peak (40 fertile females, ‘ff_40’) compared with when females were no longer fertile and were incubating (0 fertile females, ‘ff_0’). Many males were much more active during the peak fertile period.
The table below summarises activity levels when there were no fertile females.
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==0) %>%
summarise(mean_perc_active=mean(perc_active),
sd_perc_active=sd(perc_active),
mean_hours_active=24/100*mean_perc_active,
sd_hours_active=24/100*sd_perc_active) %>%
flextable()mean_perc_active | sd_perc_active | mean_hours_active | sd_hours_active |
|---|---|---|---|
34.7 | 15.6 | 8.33 | 3.74 |
The table below summarises activity levels when there were the maximum number of fertile females.
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==40) %>%
summarise(mean_perc_active=mean(perc_active),
sd_perc_active=sd(perc_active),
mean_hours_active=24/100*mean_perc_active,
sd_hours_active=24/100*sd_perc_active) %>%
flextable()mean_perc_active | sd_perc_active | mean_hours_active | sd_hours_active |
|---|---|---|---|
71.8 | 25.4 | 17.2 | 6.1 |
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
select(colour_combo, perc_active, n_FertileFemales) %>%
group_by(colour_combo, n_FertileFemales) %>%
summarise(mean_perc_active=mean(perc_active)) %>%
pivot_wider(values_from = 'mean_perc_active', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
filter(!is.na(ff_40) & !is.na(ff_0)) %>%
mutate(diff=ff_40-ff_0,
diff_hours=24/100*diff) %>%
ungroup() %>%
summarise(average_within_ind_diff_hours = mean(diff_hours),
sd_within_ind_diff_hours = sd(diff_hours),
n=n()) %>%
flextable()average_within_ind_diff_hours | sd_within_ind_diff_hours | n |
|---|---|---|
10.1 | 7.32 | 16 |
Within individual males, the average difference in activity between peak fertile females and no fertile females was over 10 hours.
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
select(colour_combo, perc_active, n_FertileFemales) %>%
group_by(n_FertileFemales) %>%
summarise(mean_perc_active=mean(perc_active)) %>%
pivot_wider(values_from = 'mean_perc_active', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
mutate(diff=ff_40-ff_0,
diff_hours=24/100*diff) %>%
ungroup() %>%
summarise(average_among_ind_diff_hours = mean(diff_hours)) %>%
flextable()average_among_ind_diff_hours |
|---|
8.91 |
Among all males, the average difference was just under 9 hours.
Prediction: Males increase their activity (% time active per day) when there are fertile females at the site.
postfertile_males <- ind_activity %>%
filter(breeding_stage=="post-fertile" & n_days>=2) %>%
group_by(bird_ID) %>%
summarise(mean(total_perc_active),
mean(mean_ODBA))
fertile_males <- ind_activity %>%
filter(breeding_stage=="fertile") %>% # already filtered to n_days>=2
group_by(bird_ID) %>%
summarise(mean(total_perc_active),
mean(mean_ODBA))Fig: Breeding stage and male % active
ind_activity_sub <- ind_activity %>%
filter(n_days>=2 & breeding_stage!="pre-fertile") %>%
mutate(breeding_stage=as.factor(breeding_stage))
ind_activity_plotdat <- rbind(ind_activity_sub %>% mutate(plot_type="box"),
ind_activity_sub %>% mutate(plot_type="point"))
ind_activity_plotdat$plot_factor <- factor(
paste(ind_activity_plotdat$breeding_stage, ind_activity_plotdat$plot_type),
levels = c("fertile box", "fertile point", "post-fertile point", "post-fertile box"))
activity_fert_boxplot <- ggplot() +
theme_classic()+
geom_boxplot(data = ind_activity_plotdat %>%
filter(plot_factor == "fertile box"),
aes(x = plot_factor,
y = total_perc_active),
outlier.shape=default_shape,
fill="white",
width=0.7) +
geom_point(data = ind_activity_plotdat %>%
filter(plot_factor == "fertile point" |
plot_factor == "post-fertile point"),
aes(x = plot_factor,
y = total_perc_active,
fill=success, colour=success, alpha=success),
size=small.point.size,
position=position_jitter(w=0.15, h=0, seed = 1)) +
geom_boxplot(data = ind_activity_plotdat %>%
filter(plot_factor == "post-fertile box"),
aes(x = plot_factor,
y = total_perc_active),
outlier.shape=default_shape,
fill="white",
width=0.7)+
geom_line(data = ind_activity_plotdat %>%
filter(plot_factor == "fertile point" |
plot_factor == "post-fertile point"),
aes(x = plot_factor,
y = total_perc_active,
colour=success, group=bird_ID, alpha=success),
size=small.point.size/3,
position=position_jitter(w=0.15, h=0, seed = 1)) +
coord_cartesian(ylim=c(0,100))+
scale_y_continuous(breaks=c(seq(0,100,by=20)),
expand=c(0,0))+
scale_x_discrete(breaks = c("fertile box", "post-fertile box"),
labels = c("fertile","post-fertile"))+
scale_colour_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
scale_fill_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
scale_alpha_manual(values=c(0.5,0.8), name="Sired offspring?")+
scale_shape_manual(values=c(default_shape), guide="none") +
xlab("\n breeding stage")+
ylab("male activity (% time)\n")+
theme(text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
legend.position="right",
legend.justification="top",
legend.background = element_blank(),
legend.box.background = element_rect(colour = "grey"))
activity_fert_boxplot
ind_activity %>%
group_by(breeding_stage, success) %>%
summarise(mean_total_perc_active = mean(total_perc_active),
sd_total_perc_active = sd(total_perc_active),
n=n(),
se_total_perc_active = sd_total_perc_active/sqrt(n)) %>%
flextable()breeding_stage | success | mean_total_perc_active | sd_total_perc_active | n | se_total_perc_active |
|---|---|---|---|---|---|
fertile | no | 58.2 | 14.9 | 65 | 1.85 |
fertile | yes | 67.1 | 12.1 | 16 | 3.03 |
post-fertile | no | 41.9 | 10.5 | 41 | 1.64 |
post-fertile | yes | 40.3 | 12.6 | 15 | 3.25 |
In the plot above, there are 81 individual males in total, PI,PI,R, O,DB,DG, PI,Y,R, Y,DB,Y, PI,PI,O, O,DB,O, DG,O,DG, DG,DG,R, DG,R,DB, R,PI,PI, Y,Y,PI, Y,Y,O, DG,O,O, Y,DG,PI, R,O,O, DB,DG,Y, O,R,PI, R,DG,Y, DG,PI,O, PI,R,PI, DB,O,DG, PI,DG,DG, Y,Y,R, Y,O,Y, Y,Y,DG, DB,Y,PI, R,DG,R, R,O,R, PI,O,R, R,R,DG, Y,O,R, DB,Y,R, R,DG,PI, O,DG,DG, DB,R,O, DG,R,DG, Y,Y,Y, DG,R,O, O,R,O, O,PI,Y, PI,R,DB, O,DB,R, DG,DB,Y, DB,DG,DG, Y,DG,DG, DB,R,R, R,DB,R, DB,PI,PI, DG,PI,Y, R,Y,Y, PI,PI,PI, DG,R,R, DB,R,Y, R,Y,R, O,Y,O, PI,DB,DB of which have post-fertile data. Males with “uncertain” success are treated as “no” success. Their inclusion or exclusion does not ever seem to affect results.
Each datapoint represents the mean activity of an individual male. The pattern we see here is similar to what was found in Lesku et al. 2012, although activity levels are shifted around 30% lower.
Model: Breeding stage and male % active
daily_activity <- daily_activity %>%
mutate(prop_activity=perc_active/100,
bird_ID=as.factor(bird_ID),
breeding_stage=as.factor(breeding_stage),
success=as.factor(success)) %>%
mutate(prop_activity=case_when(prop_activity==1 ~ 0.99999,
TRUE ~ prop_activity))stage_m <- glmmTMB(prop_activity ~ breeding_stage + (breeding_stage|bird_ID),
data=daily_activity,
#data=daily_activity %>% filter(likely_undetected_success=="no"),
family=beta_family(link="logit"))
summary(stage_m)## Family: beta ( logit )
## Formula: prop_activity ~ breeding_stage + (breeding_stage | bird_ID)
## Data: daily_activity
##
## AIC BIC logLik deviance df.resid
## -518 -490 265 -530 796
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## bird_ID (Intercept) 0.415 0.644
## breeding_stagepost-fertile 0.585 0.765 -0.93
## Number of obs: 802, groups: bird_ID, 81
##
## Dispersion parameter for beta family (): 6.78
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5391 0.0796 6.77 1.2e-11 ***
## breeding_stagepost-fertile -0.8601 0.1073 -8.02 1.1e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(stage_m, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) 0.383 0.695 0.539
## breeding_stagepost-fertile -1.070 -0.650 -0.860
## Std.Dev.(Intercept)|bird_ID 0.532 0.780 0.644
## Std.Dev.breeding_stagepost-fertile|bird_ID 0.605 0.968 0.765
## Cor.breeding_stagepost-fertile.(Intercept)|bird_ID -0.967 -0.774 -0.929simulationOutput <- simulateResiduals(stage_m)
testResiduals(simulationOutput)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.004
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.8, p-value = 0.02
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 7, observations = 802, p-value = 0.7
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
## 0.00352 0.01790
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.00873Important note We are using the DHARMA package to visualise whether models met model assumptions. This package also provides statistical tests and p-values for assumption violations, but we have focused on the visualisation of these assumptions instead. Sometimes there were statistically significant deviations from a particular assumption, but these appear reasonably minor and we believe they are of no cause for concern.
stage_m_flex <- as_flextable(stage_m)
stage_m_flexgroup | Estimate | Standard Error | statistic | p-value | |||
|---|---|---|---|---|---|---|---|
Fixed effects | |||||||
(Intercept) | 0.539 | 0.080 | 6.775 | 0.0000 | *** | ||
breeding_stagepost-fertile | -0.860 | 0.107 | -8.016 | 0.0000 | *** | ||
Random effects | |||||||
bird_ID | sd__(Intercept) | 0.644 | |||||
bird_ID | sd__breeding_stagepost-fertile | 0.765 | |||||
bird_ID | cor__(Intercept).breeding_stagepost-fertile | -0.929 | |||||
Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05 | |||||||
square root of the estimated residual variance: 6.8 | |||||||
data's log-likelihood under the model: 265.2 | |||||||
Akaike Information Criterion: -518.5 | |||||||
Bayesian Information Criterion: -490.4 | |||||||
if(save_tables==TRUE){
save_as_docx(stage_m_flex,
pr_section = sect_properties,
path = glue("{output_path}/breeding_stage_model.docx"))
}When ‘potentially’ successful males are excluded the results remain the same. There are no significant issues with model dispersion or influence of outliers.
Male spent more time active during the fertile stage, compared with the post-fertile stage. The inclusion of males that were more likely to have sired additional offspring (predated before they could be sampled) had no effect on these results.
Prediction 2 Male activity (% time active per day) increases with the number of fertile females at the site.
We also examine this trend in more detail, by looking at activity and the number of fertile females each day across the season.
Fig: Date and male % active
date_activity_n <- daily_activity %>%
group_by(local_date) %>%
summarise(n_ind=n_distinct(colour_combo))date_labels = c("",
"12 June", "", "", "",
"16 June", "", "", "",
"20 June", "", "", "",
"24 June", "", "", "",
"28 June", "", "", "",
"2 July", "","","",
"6 July", "","","",
"10 July", "")boxplot_active_date <- ggplot(daily_activity, aes(x=local_date,
y=perc_active,
group=local_date))+
# theme_classic()+
geom_boxplot(colour="gray30", outlier.shape = NA)+
geom_point(size=small.point.size,
aes(colour=success, alpha=success, fill=success),
shape=default_shape)+
geom_vline(xintercept=c(seq(min(daily_activity$local_date)-24*60*30,
max(daily_activity$local_date)+24*60*30,
by="days")),
size=0.5,
colour="darkgrey")+
geom_vline(xintercept=max(daily_activity$local_date)+24*60*30,
size=0.5,
colour="black")+
coord_flip(ylim=c(0,115),
xlim = c(max(daily_activity$local_date)+24*60*30,
min(daily_activity$local_date)-24*60*30))+
scale_x_continuous(breaks = c(seq(
local_date_to_posix('2022-06-11'),
local_date_to_posix('2022-07-11'),
by="days")),
labels=date_labels,
expand=c(0,0))+
scale_y_continuous(breaks=seq(0,100,by=20),
expand=c(0,0)) +
scale_colour_manual(values = c(fail_colour, success_colour), guide="none")+
scale_fill_manual(values = c(fail_colour, success_colour), guide="none")+
scale_alpha_manual(values = c(0.5, 0.8), guide="none")+
xlab("date\n")+
ylab("\nmale activity (% time)")+
geom_text(data = date_activity_n,
family=font.family,
size=(label.size-2)/.pt,
colour="gray50",
aes(local_date, 107, label = n_ind),
vjust = 0.3,
hjust = 0)+
theme(axis.ticks.length.y = unit(0.07, "cm"),
text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size),
axis.line.y = element_line(size=0.5, color="black"),
panel.background=element_blank())
boxplot_active_date
hor_barplot_ff <- barplot_ff +
theme(axis.line.x=element_blank()) +
coord_flip(xlim = c(max(daily_activity$local_date)+24*60*30,
min(daily_activity$local_date)-24*60*30),
ylim=c(0,42)) +
geom_vline(xintercept=c(seq(min(daily_activity$local_date)-24*60*30,
max(daily_activity$local_date)+24*60*30,
by="days")),
size=0.5,
colour="darkgrey")+
geom_vline(xintercept=max(daily_activity$local_date)+24*60*30,
size=0.5,
colour="black")+
scale_x_continuous(breaks = c(seq(
local_date_to_posix('2022-06-11'),
local_date_to_posix('2022-07-11'),
by="days")),
expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
breaks = seq(0,45,by=10))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size)) +
ylab("\n number of fertile females")
season_act_plot <- boxplot_active_date +
hor_barplot_ff +
plot_layout(
nrow=1,
widths = c(1, 0.5))
season_act_plot
Model: Male % active and number of fertile females
ff_m <- glmmTMB(prop_activity ~
scale(n_FertileFemales) +
(scale(n_FertileFemales)|bird_ID),
data = daily_activity,
family = beta_family(link="logit"))
summary(ff_m)## Family: beta ( logit )
## Formula:
## prop_activity ~ scale(n_FertileFemales) + (scale(n_FertileFemales) |
## bird_ID)
## Data: daily_activity
##
## AIC BIC logLik deviance df.resid
## -590 -562 301 -602 796
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## bird_ID (Intercept) 0.186 0.432
## scale(n_FertileFemales) 0.176 0.419 0.77
## Number of obs: 802, groups: bird_ID, 81
##
## Dispersion parameter for beta family (): 7.61
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2660 0.0562 4.73 0.000002228925850 ***
## scale(n_FertileFemales) 0.4326 0.0564 7.67 0.000000000000018 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(ff_m, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) 0.156 0.376 0.266
## scale(n_FertileFemales) 0.322 0.543 0.433
## Std.Dev.(Intercept)|bird_ID 0.350 0.533 0.432
## Std.Dev.scale(n_FertileFemales)|bird_ID 0.335 0.524 0.419
## Cor.scale(n_FertileFemales).(Intercept)|bird_ID 0.471 0.886 0.774simulationOutput <- simulateResiduals(ff_m)
testResiduals(simulationOutput)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.01
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9, p-value = 0.04
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
## 0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.0125ff_flex <- as_flextable(ff_m)
if(save_tables==TRUE){
save_as_docx(ff_flex,
pr_section = sect_properties,
path = glue("{output_path}/fertile_females_model.docx"))
}Male activity (% time) was higher when there were more fertile females at the site.
We can also look at a variation of this model to better understand within- and among-individual variation.
Model: Within- and among-individual variation in male activity
mean_fertile_females <- daily_activity %>%
group_by(colour_combo, bird_ID) %>%
summarise(ind_mean_ff = mean(n_FertileFemales))
daily_activity <- left_join(daily_activity, mean_fertile_females,
join_by('colour_combo','bird_ID')) %>%
mutate(ind_dev_ff=n_FertileFemales-ind_mean_ff)In our study, there are two proximate mechanisms affecting the relationship between number of fertile females and male activity: individual changes in activity (within-individual variation) and differences in the presence of males at the site (between-individual variation). One method for distinguishing between these two mechanisms is within-subject centering (van de Pol and Wright 2009; Dingemanse and Dochtermann 2013).
From Dingemanse and Dochtermann 2013: “[The] conflation of between- and within-individual effects can be separated by calculating the mean covariate value for the individual, and the deviation of the covariate from the mean for an individual for each measurement. Both of these predictor variables are then included in modelling the phenotypic response”.
The effect of the latter (the “deviation” term) tells us the effect of the number of fertile females on male activity. The former (the “individual mean” term) tells us the effect of individuals being measured when there are different numbers of fertile females.
We can again include a random slope, so that we can also look at 1) among-individual variation in the response to changes in the number of fertile females, and 2) how activity ‘elevation’ relates to activity ‘slope’ (covariance between random intercept and slope).
ff_m2 <- glmmTMB(prop_activity ~ ind_mean_ff + ind_dev_ff +
(ind_dev_ff|bird_ID),
data = daily_activity,
family = beta_family(link="logit"))
summary(ff_m2)## Family: beta ( logit )
## Formula:
## prop_activity ~ ind_mean_ff + ind_dev_ff + (ind_dev_ff | bird_ID)
## Data: daily_activity
##
## AIC BIC logLik deviance df.resid
## -591 -558 302 -605 795
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## bird_ID (Intercept) 0.246 0.4956
## ind_dev_ff 0.001 0.0317 0.79
## Number of obs: 802, groups: bird_ID, 81
##
## Dispersion parameter for beta family (): 7.71
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.21055 0.20567 1.02 0.31
## ind_mean_ff 0.00417 0.00786 0.53 0.60
## ind_dev_ff 0.03480 0.00446 7.80 6e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(ff_m2, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -0.1926 0.6137 0.21055
## ind_mean_ff -0.0112 0.0196 0.00417
## ind_dev_ff 0.0261 0.0435 0.03480
## Std.Dev.(Intercept)|bird_ID 0.4054 0.6060 0.49564
## Std.Dev.ind_dev_ff|bird_ID 0.0249 0.0403 0.03168
## Cor.ind_dev_ff.(Intercept)|bird_ID 0.5369 0.8885 0.78956simulationOutput <- simulateResiduals(ff_m)
testResiduals(simulationOutput)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.01
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9, p-value = 0.04
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
## 0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.0125Based on this model, we can make the following conclusions:
While both within- and between-individual variation contribute to the relationship between mating opportunities and male activity, within-individual variation has a stronger effect
There does not seem to be that much among-individual variation in activity slopes (i.e. most males are increasing their activity by a similar amount as number of fertile females increases?)
There is a fairly strong positive correlation between individual ‘elevation’ and ‘slope’. So, males who were more active at their ‘centred’ number of fertile females are changing their activity more than other males.
ff2_flex <- as_flextable(ff_m2)
if(save_tables==TRUE){
save_as_docx(ff2_flex,
pr_section = sect_properties,
path = glue("{output_path}/fertile_females_variation_model.docx"))
}Fig: Male % active and number of fertile females
m_effects = effects::allEffects(ff_m, xlevels = 160)
m_ff_effects = m_effects$`scale(n_FertileFemales)` %>% data.frameactivity_ff_plot <- ggplot(daily_activity, aes(x=n_FertileFemales,
y=perc_active))+
theme_classic()+
#geom_boxplot(colour="gray30", outlier.shape = NA)+
geom_point(size=small.point.size,
shape=default_shape,
alpha=0.15,
colour=male_colour,
fill=male_colour)+
geom_line(data = m_ff_effects,
aes(y = fit*100, x = n_FertileFemales),
colour=col_scale[1]) +
geom_ribbon(data = m_ff_effects,
aes(y = fit*100, x = n_FertileFemales, ymin = lower*100, ymax = upper*100),
alpha = 0.2, fill=col_scale[1])+
coord_cartesian(ylim = c(0,102)) +
scale_y_continuous(breaks=seq(0,100,by=20),
expand = c(0,0)) +
xlab("\nnumber of fertile females\n")+
ylab("daily activity (% time)\n")+
theme(axis.ticks.length.y = unit(0.07, "cm"),
text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
activity_ff_plot
if(save_plots==TRUE){
ggsave(plot=activity_ff_plot,
filename='figure-activity-fertility-means-plot.jpg',
path=output_path,
width=16,
height=12.5,
units="cm",
dpi=300)
}Combined figure: Date and male activity
layout <- "
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
"
fig1 <- activity_fert_boxplot + ylab("\nmale activity (% time)\n") + xlab("\nbreeding stage\n\n") +
plot_spacer() +
#plot_spacer() +
season_act_plot +
plot_layout(design=layout) +
plot_layout(nrow=2)+
plot_annotation(tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')') &
theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"))
fig1
if(save_plots==TRUE){
ggsave(plot=fig1,
filename='figure1-activity-stage_date.jpg',
path=output_path,
width=16,
height=20,
units="cm",
dpi=300)
}Q: Do males increase the intensity of their activity when there are more mating opportunities?
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
select(colour_combo, mean_ODBA, n_FertileFemales) %>%
group_by(colour_combo, n_FertileFemales) %>%
summarise(mean_ODBA=mean(mean_ODBA)) %>%
pivot_wider(values_from = 'mean_ODBA', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
filter(!is.na(ff_40) & !is.na(ff_0)) %>%
mutate(diff=ff_40-ff_0) %>%
flextable()colour_combo | ff_40 | ff_0 | diff |
|---|---|---|---|
DB,O,DG | 0.422 | 0.201 | 0.22086 |
DB,PI,PI | 0.193 | 0.189 | 0.00373 |
DB,R,O | 0.566 | 0.138 | 0.42849 |
DG,R,DG | 0.211 | 0.188 | 0.02260 |
O,DB,DG | 0.530 | 0.211 | 0.31961 |
O,DG,DG | 0.505 | 0.180 | 0.32515 |
O,PI,Y | 0.164 | 0.182 | -0.01871 |
PI,DG,DG | 0.596 | 0.205 | 0.39034 |
PI,PI,PI | 0.552 | 0.160 | 0.39183 |
PI,Y,R | 0.193 | 0.175 | 0.01812 |
R,DB,R | 0.390 | 0.240 | 0.15035 |
R,Y,R | 0.315 | 0.117 | 0.19782 |
Y,DB,Y | 0.449 | 0.193 | 0.25670 |
Y,DG,PI | 0.248 | 0.177 | 0.07112 |
Y,O,R | 0.501 | 0.175 | 0.32574 |
Y,Y,DG | 0.498 | 0.204 | 0.29401 |
The table above shows the mean ODBA of males (each with a unique colour_combo) when female fertility was at its peak (40 fertile females, ‘ff_40’) compared with when females were no longer fertile and were incubating (0 fertile females, ‘ff_0’). Almost all males showed more intense activity during the peak fertile period.
The table below summarises mean ODBA levels when there were no fertile females.
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==0) %>%
summarise(sd_ODBA=sd(mean_ODBA),
mean_ODBA=mean(mean_ODBA)) %>%
flextable()sd_ODBA | mean_ODBA |
|---|---|
0.0343 | 0.186 |
The table below summarises mean ODBA levels when there were the maximum number of fertile females.
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales==40) %>%
summarise(sd_ODBA=sd(mean_ODBA),
mean_ODBA=mean(mean_ODBA)) %>%
flextable()sd_ODBA | mean_ODBA |
|---|---|
0.168 | 0.388 |
Prediction: Males increase the intensity of their activity (mean ODBA) when there are fertile females at the site.
Fig: Breeding stage and male ODBA
odba_fert_boxplot <- ggplot() +
theme_classic()+
geom_boxplot(data = ind_activity_plotdat %>%
filter(plot_factor == "fertile box"),
aes(x = plot_factor,
y = mean_ODBA),
outlier.shape=default_shape,
fill="white",
width=0.7) +
geom_point(data = ind_activity_plotdat %>%
filter(plot_factor == "fertile point" |
plot_factor == "post-fertile point"),
aes(x = plot_factor,
y = mean_ODBA,
fill=success, colour=success, alpha=success),
size=small.point.size,
position=position_jitter(w=0.15, h=0, seed = 1)) +
geom_boxplot(data = ind_activity_plotdat %>%
filter(plot_factor == "post-fertile box"),
aes(x = plot_factor,
y = mean_ODBA),
outlier.shape=default_shape,
fill="white",
width=0.7)+
geom_line(data = ind_activity_plotdat %>%
filter(plot_factor == "fertile point" |
plot_factor == "post-fertile point"),
aes(x = plot_factor,
y = mean_ODBA,
colour=success, group=bird_ID, alpha=success),
size=small.point.size/3,
position=position_jitter(w=0.15, h=0, seed = 1)) +
coord_cartesian(ylim=c(0,0.6))+
scale_y_continuous(breaks=c(seq(0,0.6,by=0.1)),
expand=c(0,0))+
scale_x_discrete(breaks = c("fertile box", "post-fertile box"),
labels = c("fertile","post-fertile"))+
scale_colour_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
scale_fill_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
scale_alpha_manual(values=c(0.5,0.8), name="Sired offspring?")+
scale_shape_manual(values=c(default_shape), guide="none") +
xlab("\n breeding stage")+
ylab("mean male ODBA (g)\n")+
theme(text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
legend.position="right",
legend.justification="top",
legend.background = element_blank(),
legend.box.background = element_rect(colour = "grey"))
odba_fert_boxplot
ind_activity %>%
group_by(breeding_stage, success) %>%
summarise(mean_mean_ODBA = mean(mean_ODBA),
sd_ODBA = sd(mean_ODBA),
n=n(),
se_ODBA = sd_ODBA/sqrt(n)) %>%
flextable()breeding_stage | success | mean_mean_ODBA | sd_ODBA | n | se_ODBA |
|---|---|---|---|---|---|
fertile | no | 0.284 | 0.0964 | 65 | 0.01196 |
fertile | yes | 0.354 | 0.0725 | 16 | 0.01811 |
post-fertile | no | 0.201 | 0.0347 | 41 | 0.00542 |
post-fertile | yes | 0.222 | 0.0484 | 15 | 0.01248 |
In the plot above, there are 81 individual males in total, PI,PI,R, O,DB,DG, PI,Y,R, Y,DB,Y, PI,PI,O, O,DB,O, DG,O,DG, DG,DG,R, DG,R,DB, R,PI,PI, Y,Y,PI, Y,Y,O, DG,O,O, Y,DG,PI, R,O,O, DB,DG,Y, O,R,PI, R,DG,Y, DG,PI,O, PI,R,PI, DB,O,DG, PI,DG,DG, Y,Y,R, Y,O,Y, Y,Y,DG, DB,Y,PI, R,DG,R, R,O,R, PI,O,R, R,R,DG, Y,O,R, DB,Y,R, R,DG,PI, O,DG,DG, DB,R,O, DG,R,DG, Y,Y,Y, DG,R,O, O,R,O, O,PI,Y, PI,R,DB, O,DB,R, DG,DB,Y, DB,DG,DG, Y,DG,DG, DB,R,R, R,DB,R, DB,PI,PI, DG,PI,Y, R,Y,Y, PI,PI,PI, DG,R,R, DB,R,Y, R,Y,R, O,Y,O, PI,DB,DB of which have post-fertile data. Males with “uncertain” success are treated as “no” success.
Each datapoint represents the mean ODBA of an individual male.
Model: Breeding stage and male ODBA
odba_stage_m <- glmmTMB(mean_ODBA ~ breeding_stage + (breeding_stage|bird_ID),
data=daily_activity,
#data=daily_activity %>% filter(likely_undetected_success=="no"),
family=beta_family(link="logit"))
summary(odba_stage_m)## Family: beta ( logit )
## Formula: mean_ODBA ~ breeding_stage + (breeding_stage | bird_ID)
## Data: daily_activity
##
## AIC BIC logLik deviance df.resid
## -1690 -1662 851 -1702 796
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## bird_ID (Intercept) 0.185 0.430
## breeding_stagepost-fertile 0.111 0.333 -0.99
## Number of obs: 802, groups: bird_ID, 81
##
## Dispersion parameter for beta family (): 32.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8712 0.0511 -17.03 <2e-16 ***
## breeding_stagepost-fertile -0.4526 0.0493 -9.18 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(odba_stage_m, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -0.971 -0.771 -0.871
## breeding_stagepost-fertile -0.549 -0.356 -0.453
## Std.Dev.(Intercept)|bird_ID 0.361 0.512 0.430
## Std.Dev.breeding_stagepost-fertile|bird_ID 0.258 0.431 0.333
## Cor.breeding_stagepost-fertile.(Intercept)|bird_ID -0.999 0.991 -0.986simulationOutput <- simulateResiduals(odba_stage_m)
testResiduals(simulationOutput)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.05, p-value = 0.02
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 6, observations = 802, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
## 0.00275 0.01621
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.00748odba_stage_m_flex <- as_flextable(odba_stage_m)
odba_stage_m_flexgroup | Estimate | Standard Error | statistic | p-value | |||
|---|---|---|---|---|---|---|---|
Fixed effects | |||||||
(Intercept) | -0.871 | 0.051 | -17.034 | 0.0000 | *** | ||
breeding_stagepost-fertile | -0.453 | 0.049 | -9.175 | 0.0000 | *** | ||
Random effects | |||||||
bird_ID | sd__(Intercept) | 0.430 | |||||
bird_ID | sd__breeding_stagepost-fertile | 0.333 | |||||
bird_ID | cor__(Intercept).breeding_stagepost-fertile | -0.986 | |||||
Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05 | |||||||
square root of the estimated residual variance: 32.4 | |||||||
data's log-likelihood under the model: 850.8 | |||||||
Akaike Information Criterion: -1,689.7 | |||||||
Bayesian Information Criterion: -1,661.6 | |||||||
if(save_tables==TRUE){
save_as_docx(odba_stage_m_flex,
pr_section = sect_properties,
path = glue("{output_path}/breeding_stage_odba_model.docx"))
}When ‘potentially’ successful males are excluded the conclusions remain the same.
as.data.frame(emmeans(odba_stage_m, ~ breeding_stage))## breeding_stage emmean SE df asymp.LCL asymp.UCL
## fertile -0.871 0.0511 Inf -0.971 -0.771
## post-fertile -1.324 0.0307 Inf -1.384 -1.264
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95odba_stage_emm_flex <- as_flextable(as.data.frame(emmeans(odba_stage_m, ~ breeding_stage)))
if(save_tables==TRUE){
save_as_docx(odba_stage_emm_flex,
pr_section = sect_properties,
path = glue("{output_path}/breeding_stage_emm.docx"))
}Male tended to show more intense activity during the fertile stage, compared with the post-fertile stage. The inclusion of males that were more likely to have sired additional offspring (predated before they could be sampled) had no effect on these results.
Prediction 2 Intensity of male activity (mean ODBA) increases with the number of fertile females at the site.
We also examine this trend in more detail, by looking at activity and the number of fertile females each day across the season.
Fig: Date and male ODBA
boxplot_odba_date <- ggplot(daily_activity, aes(x=local_date,
y=mean_ODBA,
group=local_date))+
geom_boxplot(colour="gray30", outlier.shape = NA)+
geom_point(size=small.point.size,
aes(colour=success, alpha=success, fill=success),
shape=default_shape)+
geom_vline(xintercept=c(seq(min(daily_activity$local_date)-24*60*30,
max(daily_activity$local_date)+24*60*30,
by="days")),
size=0.5,
colour="darkgrey")+
geom_vline(xintercept=max(daily_activity$local_date)+24*60*30,
size=0.5,
colour="black")+
coord_flip(ylim=c(0,0.9),
xlim = c(max(daily_activity$local_date)+24*60*30,
min(daily_activity$local_date)-24*60*30))+
scale_x_continuous(breaks = c(seq(
local_date_to_posix('2022-06-11'),
local_date_to_posix('2022-07-11'),
by="days")),
labels=date_labels,
expand=c(0,0))+
scale_y_continuous(breaks=seq(0,0.9,by=0.2),
expand=c(0,0)) +
scale_colour_manual(values = c(fail_colour, success_colour), guide="none")+
scale_fill_manual(values = c(fail_colour, success_colour), guide="none")+
scale_alpha_manual(values = c(0.5, 0.8), guide="none")+
xlab("date\n")+
ylab("\nmean male ODBA (g)")+
geom_text(data = date_activity_n,
family=font.family,
size=(label.size-2)/.pt,
colour="gray50",
aes(local_date, 107, label = n_ind),
vjust = 0.3,
hjust = 0)+
theme(axis.ticks.length.y = unit(0.07, "cm"),
text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size),
axis.line.y = element_line(size=0.5, color="black"),
panel.background=element_blank())
boxplot_odba_date
season_odba_plot <- boxplot_odba_date +
hor_barplot_ff +
plot_layout(
nrow=1,
widths = c(1, 0.5))
season_odba_plot
Model: Fertile females and ODBA
ff_m_odba <- glmmTMB(mean_ODBA ~ scale(n_FertileFemales) +
(scale(n_FertileFemales)|bird_ID),
data = daily_activity,
family = beta_family(link="logit"))
summary(ff_m_odba)## Family: beta ( logit )
## Formula:
## mean_ODBA ~ scale(n_FertileFemales) + (scale(n_FertileFemales) | bird_ID)
## Data: daily_activity
##
## AIC BIC logLik deviance df.resid
## -1764 -1736 888 -1776 796
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## bird_ID (Intercept) 0.1269 0.356
## scale(n_FertileFemales) 0.0412 0.203 0.63
## Number of obs: 802, groups: bird_ID, 81
##
## Dispersion parameter for beta family (): 39.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0219 0.0429 -23.80 <2e-16 ***
## scale(n_FertileFemales) 0.2398 0.0287 8.36 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(ff_m_odba, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -1.106 -0.938 -1.022
## scale(n_FertileFemales) 0.184 0.296 0.240
## Std.Dev.(Intercept)|bird_ID 0.298 0.426 0.356
## Std.Dev.scale(n_FertileFemales)|bird_ID 0.159 0.259 0.203
## Cor.scale(n_FertileFemales).(Intercept)|bird_ID 0.318 0.788 0.628simulationOutput <- simulateResiduals(ff_m_odba)
testResiduals(simulationOutput)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.008
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.2
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
## 0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.0125ff_odba_flex <- as_flextable(ff_m_odba)
if(save_tables==TRUE){
save_as_docx(ff_odba_flex,
pr_section = sect_properties,
path = glue("{output_path}/fertile_females_odba_model.docx"))
}Model: Within- and among-individual variation
ff_m_odba2 <- glmmTMB(mean_ODBA ~ ind_mean_ff + ind_dev_ff +
(ind_dev_ff|bird_ID),
data = daily_activity,
family = beta_family(link="logit"))
summary(ff_m_odba2)## Family: beta ( logit )
## Formula: mean_ODBA ~ ind_mean_ff + ind_dev_ff + (ind_dev_ff | bird_ID)
## Data: daily_activity
##
## AIC BIC logLik deviance df.resid
## -1768 -1736 891 -1782 795
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## bird_ID (Intercept) 0.143905 0.379
## ind_dev_ff 0.000226 0.015 0.68
## Number of obs: 802, groups: bird_ID, 81
##
## Dispersion parameter for beta family (): 39.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03498 0.15620 -6.63 0.000000000034 ***
## ind_mean_ff 0.00197 0.00592 0.33 0.74
## ind_dev_ff 0.01817 0.00218 8.34 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(ff_m_odba2, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -1.34113 -0.7288 -1.03498
## ind_mean_ff -0.00964 0.0136 0.00197
## ind_dev_ff 0.01391 0.0224 0.01817
## Std.Dev.(Intercept)|bird_ID 0.31889 0.4513 0.37935
## Std.Dev.ind_dev_ff|bird_ID 0.01172 0.0193 0.01504
## Cor.ind_dev_ff.(Intercept)|bird_ID 0.41008 0.8132 0.67846simulationOutput <- simulateResiduals(ff_m_odba2)
testResiduals(simulationOutput)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.05, p-value = 0.04
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.4
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
## 0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.0125Based on this model, we can make the following conclusions:
The relationship between mating opportunities and male activity is explained by within-individual variation
There does not seem to be that much among-individual variation in activity slopes (i.e. most males are increasing their activity by a similar amount as number of fertile females increases?)
There is a fairly strong negative correlation between individual ‘elevation’ and ‘slope’, like what we saw with percent time active
ff_odba2_flex <- as_flextable(ff_m_odba2)
if(save_tables==TRUE){
save_as_docx(ff_odba2_flex,
pr_section = sect_properties,
path = glue("{output_path}/fertile_females_odba_variation_model.docx"))
}Fig: Male ODBA and number of fertile females
m_effects = effects::allEffects(ff_m_odba, xlevels = 160)
m_ff_odba_effects = m_effects$`scale(n_FertileFemales)` %>% data.frame
odba_ff_plot <- ggplot(daily_activity, aes(x=n_FertileFemales,
y=mean_ODBA))+
theme_classic()+
#geom_boxplot(colour="gray30", outlier.shape = NA)+
geom_point(size=small.point.size,
shape=default_shape,
alpha=0.15,
colour=male_colour,
fill=male_colour)+
geom_line(data = m_ff_odba_effects,
aes(y = fit, x = n_FertileFemales),
colour=col_scale[1]) +
geom_ribbon(data = m_ff_odba_effects,
aes(y = fit, x = n_FertileFemales, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[1])+
coord_cartesian(ylim = c(0,0.9)) +
scale_y_continuous(breaks=seq(0,0.9,by=0.2),
expand = c(0,0)) +
xlab("\nnumber of fertile females\n")+
ylab("mean daily ODBA (g)\n")+
theme(axis.ticks.length.y = unit(0.07, "cm"),
text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
odba_ff_plot
daily_activity %>%
ungroup() %>%
filter(n_FertileFemales == 0 | n_FertileFemales==40) %>%
group_by(n_FertileFemales) %>%
summarise(mean_mean_ODBA = mean(mean_ODBA),
sd_odba_active = sd(mean_ODBA),
n=n(),
se_odba_active=sd_odba_active/sqrt(n))## # A tibble: 2 × 5
## n_FertileFemales mean_mean_ODBA sd_odba_active n se_odba_active
## <int> <dbl> <dbl> <int> <dbl>
## 1 0 0.186 0.0343 73 0.00402
## 2 40 0.388 0.168 122 0.0152Combined figure: Date and male ODBA
layout <- "
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
"
fig_s1 <- odba_fert_boxplot + ylab("\nmean male ODBA (g)\n") + xlab("\nbreeding stage\n\n") +
plot_spacer() +
#plot_spacer() +
season_odba_plot +
plot_layout(design=layout) +
plot_layout(nrow=2)+
plot_annotation(tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')') &
theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"))
fig_s1
if(save_plots==TRUE){
ggsave(plot=fig_s1,
filename='figure-odba-stage_date.jpg',
path=output_path,
width=16,
height=20,
units="cm",
dpi=300)
}Combined figure: Fertile females and male activity
fig2 <- activity_ff_plot +
plot_spacer() +
odba_ff_plot +
plot_layout(nrow = 1,
widths=c(1,0.1,1),
heights=c(1))+
plot_annotation(tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')') &
theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
plot.tag.position = c(.27, .96))
fig2
if(save_plots==TRUE){
ggsave(plot=fig2,
filename='figure2-activity-fertilefemales.jpg',
path=output_path,
width=16,
height=7,
units="cm",
dpi=300)
}rm(stage_m, stage_m_flex,
m_ff_odba_effects, m_ff_effects, m_effects,
ff_flex, ff_m, ff_m_odba, ff_m2,
ff_odba_flex, ff2_flex)Q: How does male behaviour relate to mating and reproductive success?
Prediction Increased time spent active (%) is associated with increased likelihood of mating success Prediction Males who spend more days at the site during the fertile period have a higher probability of mating success Prediction Males with higher mean ODBA have a higher probability of mating success.
In the plots and analyses below, we only include data from the fertile stage, and only for males with at least 2 days of data during the fertile stage.
However, it’s a problem to exclude males with limited data for the tenure analysis, because then we’re excluding all the males who stayed for less than 2. So we analyse the effects of tenure on mating success in a separate model, including the full dataset.
ind_activity <- ind_activity %>%
mutate(prop_activity=total_perc_active/100,
recording_hours=n_windows*window_length/60,
arrival_day=as.numeric(difftime(arrival_date,
local_date_to_posix("2022-06-01"),
units="days")),
success_num = case_when(success=="yes" ~ 1,
success=="no" ~ 0))
ind_fertile_activity <- ind_activity %>%
filter(breeding_stage=="fertile" & n_days>=min_acc_days)
ind_tenure <- left_join(male_tenure, father_summary,
join_by('colour_combo'=='LL.father',
'bird_ID'=='ID.father')) %>%
mutate(arrival_day=as.numeric(difftime(arrival_date,
local_date_to_posix("2022-06-01"),
units="days")),
likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
TRUE ~ "less likely"),) %>%
replace(is.na(.), 0)ind_fertile_activity %>%
ungroup() %>%
summarise(min_total_perc_active=min(total_perc_active),
max_total_perc_active=max(total_perc_active),
mean_total_perc_active=mean(total_perc_active),
sd_total_perc_active=sd(total_perc_active),
n=n(),
se_total_perc_active=sd_total_perc_active/sqrt(n)) %>%
flextable()min_total_perc_active | max_total_perc_active | mean_total_perc_active | sd_total_perc_active | n | se_total_perc_active |
|---|---|---|---|---|---|
15.4 | 84.5 | 59.9 | 14.8 | 81 | 1.64 |
The values below show the correlation between variables considered for models.
cor(ind_fertile_activity$tenure_days, ind_fertile_activity$recording_hours)## [1] 0.749cor(ind_fertile_activity$tenure_days, ind_fertile_activity$prop_activity)## [1] 0.19cor(ind_fertile_activity$prop_activity, ind_fertile_activity$arrival_day)## [1] -0.228cor(ind_fertile_activity$tenure_days, ind_fertile_activity$arrival_day)## [1] -0.201Tenure is highly correlated with number of hours of data recorded, which makes sense. So we do not include both these variables in the model. Tenure is the more interesting and biologically meaningful variable, so we focus on that.
Also note that mean ODBA and proportion of time spent active are correlated.
cor.test(ind_fertile_activity$mean_ODBA, ind_fertile_activity$prop_activity)##
## Pearson's product-moment correlation
##
## data: ind_fertile_activity$mean_ODBA and ind_fertile_activity$prop_activity
## t = 19, df = 79, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.853 0.937
## sample estimates:
## cor
## 0.903Mating success (number of females with which a male sired offspring)
Model: % active and mating success
There are often concerns about models with Poisson distributions. We confirmed that all model assumptions were met and that there were no issues with overdispersion. Running the model with a Gaussian distribution violates model assumptions but still produces the same conclusions.
success_m1 <- glmmTMB(n_nests ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day),
data=ind_fertile_activity,
# data=ind_fertile_activity %>% filter(likely_undetected_success=="no"),
family=poisson)
summary(success_m1)## Family: poisson ( log )
## Formula:
## n_nests ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
##
## AIC BIC logLik deviance df.resid
## 97.8 107.4 -44.9 89.8 77
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.694 0.293 -5.78 0.0000000073 ***
## scale(prop_activity) 0.436 0.264 1.66 0.098 .
## scale(log(tenure_days)) 0.554 0.332 1.67 0.095 .
## scale(arrival_day) 0.109 0.253 0.43 0.666
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(success_m1, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -2.2682 -1.120 -1.694
## scale(prop_activity) -0.0802 0.953 0.436
## scale(log(tenure_days)) -0.0956 1.204 0.554
## scale(arrival_day) -0.3874 0.606 0.109simulationOutput <- simulateResiduals(fittedModel = success_m1, plot=F)
testResiduals(simulationOutput)
## $uniformity
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.07, p-value = 0.8
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.7
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 81, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0000 0.0123
## sample estimates:
## outlier frequency (expected: 0.00197530864197531 )
## 0success_m1_flex <- as_flextable(success_m1)
if(save_tables==TRUE){
save_as_docx(success_m1_flex,
pr_section = sect_properties,
path = glue("{output_path}/activity_mating_model.docx"))
}The effect of activity is similar when the males with territories around predated (unsampled) nests are excluded, although effect of tenure becomes larger.
m1_effects = effects::allEffects(success_m1, xlevels = 100)
m1_active_eff = m1_effects$`scale(prop_activity)` %>% data.frame
m1_tenure_eff = m1_effects$`scale(log(tenure_days))` %>% data.frame
m1_effects2 = effects::allEffects(success_m1, xlevels = 100, confidence.level=0.95)
m1_active_eff2 = m1_effects2$`scale(prop_activity)` %>% data.framem1_effects = as.data.frame(confint(success_m1, level=0.95))
write.csv(m1_effects,
file=glue("{output_supp_effects}/activity-mating-effects_{window_length}min_{data_type}-{threshold_type}thresh.csv"))Fig: % active and mating success
success_active_plot <- ggplot(ind_fertile_activity)+
geom_line(data = m1_active_eff,
aes(y = fit, x = prop_activity*100),
colour=col_scale[2]) +
geom_ribbon(data = m1_active_eff,
aes(y = fit, x =prop_activity*100, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[2])+
geom_jitter(aes(x=prop_activity*100, y=n_nests, shape=likely_extra_nests),
fill=male_colour,
colour=male_colour,
size=point_size,
alpha=point_alpha,
width=0.005, height=0.04) +
theme_classic()+
scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,100,by=20))) +
scale_y_continuous(breaks=c(seq(0,2,by=1))) +
coord_cartesian(xlim=c(10,90), ylim=c(0,2.1))+
xlab("\n activity (% time)") +
ylab("mating success\n")+
theme(legend.position="right",
plot.margin=unit(c(0,0,0,0), 'lines'),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
success_active_plot
Model: Tenure and mating success
success_mt <- glmmTMB(n_nests ~ scale(log(tenure_days_exact)) + scale(arrival_day),
data=ind_tenure,
# data=ind_tenure %>% filter(potential_n_extra_nests==0),
family=poisson)
summary(success_mt)## Family: poisson ( log )
## Formula: n_nests ~ scale(log(tenure_days_exact)) + scale(arrival_day)
## Data: ind_tenure
##
## AIC BIC logLik deviance df.resid
## 99.7 107.5 -46.9 93.7 96
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1857 0.3982 -5.49 0.00000004 ***
## scale(log(tenure_days_exact)) 1.4240 0.5590 2.55 0.011 *
## scale(arrival_day) 0.0242 0.2527 0.10 0.924
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(success_mt, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -2.966 -1.405 -2.1857
## scale(log(tenure_days_exact)) 0.328 2.520 1.4240
## scale(arrival_day) -0.471 0.519 0.0242simulationOutput <- simulateResiduals(fittedModel = success_mt, plot=F)
plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.7
## alternative hypothesis: two.sidedsuccess_mt_flex <- as_flextable(success_mt)
if(save_tables==TRUE){
save_as_docx(success_mt_flex,
pr_section = sect_properties,
path = glue("{output_path}/tenure_mating_model.docx"))
}The effect of activity is similar when the males with territories around predated (unsampled) nests are excluded, although effect of tenure becomes larger.
mt_effects = effects::allEffects(success_mt, xlevels = 100)
mt_tenure_eff = mt_effects$`scale(log(tenure_days_exact))` %>% data.frameFig: Tenure and mating success
success_tenure_plot <- ggplot(ind_tenure)+
geom_line(data = mt_tenure_eff,
aes(y = fit, x = tenure_days_exact),
colour=col_scale[1]) +
geom_ribbon(data = mt_tenure_eff,
aes(y = fit, x =tenure_days_exact, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[1]) +
geom_jitter(aes(x=tenure_days_exact, y=n_nests, shape=likely_extra_nests),
fill=male_colour,
colour=male_colour,
size=point_size,
alpha=point_alpha,
width=0.01, height=0.05) +
theme_classic()+
scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,30,by=10))) +
coord_cartesian(xlim=c(0,32), ylim=c(0,2.1))+
scale_y_continuous(breaks=c(seq(0,2,by=1))) +
xlab("\n tenure (days)") +
ylab("mating success\n")+
theme(legend.position="right",
plot.margin=unit(c(0,0,0,0), 'lines'),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
success_tenure_plot
I double-checked these long tenures of unsuccessful males - they were mostly at the site or very close by, including during the fertile period. In other words, they had plenty of opportunity to sire offspring at the study site, but didn’t.
ind_tenure %>%
ungroup() %>%
summarise(mean_tenure=mean(tenure_days_exact),
sd_tenure=sd(tenure_days_exact),
n=n(),
se_tenure=sd_tenure/sqrt(n))## mean_tenure sd_tenure n se_tenure
## 1 9.49 6.25 99 0.628ind_fertile_activity %>%
filter(success=="yes") %>%
select(tenure_days_exact, n_nests) %>%
arrange(tenure_days_exact)## # A tibble: 16 × 10
## # Groups: bird_ID, colour_combo, breeding_stage, body_mass, n_eggs, n_nests,
## # success, likely_undetected_success, potential_n_extra_nests [16]
## bird_ID colour_combo breeding_stage body_mass n_eggs success
## <int> <chr> <chr> <dbl> <int> <chr>
## 1 141269127 DG,DB,DG fertile 96.6 4 yes
## 2 141269182 O,Y,O fertile 97.9 2 yes
## 3 141269175 DG,R,R fertile 98.5 3 yes
## 4 141269043 DG,O,O fertile 101. 4 yes
## 5 141269152 Y,Y,Y fertile 100. 1 yes
## 6 141269146 R,DG,PI fertile 102. 2 yes
## 7 141269046 DB,DG,Y fertile 95.3 4 yes
## 8 141269169 DB,R,R fertile 97.5 4 yes
## 9 141269026 O,DB,O fertile 94.7 4 yes
## 10 141269125 Y,Y,R fertile 95.8 1 yes
## 11 141269144 DB,Y,R fertile 102. 3 yes
## 12 141269147 O,DG,DG fertile 105. 5 yes
## 13 141269148 DB,R,O fertile 108. 3 yes
## 14 141269044 Y,DG,PI fertile 92.5 5 yes
## 15 141269114 PI,R,PI fertile 106 3 yes
## 16 141269009 PI,PI,R fertile 101. 6 yes
## # ℹ 4 more variables: likely_undetected_success <chr>,
## # potential_n_extra_nests <int>, tenure_days_exact <dbl>, n_nests <int>Model: ODBA and mating success
success_m2 <- glmmTMB(n_nests ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day),
data=ind_fertile_activity,
# data=ind_fertile_activity %>% filter(likely_undetected_success=="no"),
family=poisson)
#family=gaussian)
summary(success_m2)## Family: poisson ( log )
## Formula:
## n_nests ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
##
## AIC BIC logLik deviance df.resid
## 95.2 104.8 -43.6 87.2 77
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.786 0.319 -5.59 0.000000023 ***
## scale(mean_ODBA) 0.624 0.283 2.21 0.027 *
## scale(log(tenure_days)) 0.585 0.353 1.66 0.098 .
## scale(arrival_day) 0.257 0.269 0.96 0.338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(success_m2, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -2.4121 -1.160 -1.786
## scale(mean_ODBA) 0.0697 1.179 0.624
## scale(log(tenure_days)) -0.1072 1.277 0.585
## scale(arrival_day) -0.2692 0.784 0.257simulationOutput <- simulateResiduals(fittedModel = success_m2, plot=F)
testResiduals(simulationOutput)
## $uniformity
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.07, p-value = 0.8
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.8
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 81, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0000 0.0123
## sample estimates:
## outlier frequency (expected: 0.00172839506172839 )
## 0plot(simulationOutput)
success_m2_flex <- as_flextable(success_m2)
if(save_tables==TRUE){
save_as_docx(success_m2_flex,
pr_section = sect_properties,
path = glue("{output_path}/odba_mating_model.docx"))
}Similar effects are observed when males near predated nests are excluded.
m2_effects = effects::allEffects(success_m2, xlevels = 100)
m2_active_eff = m2_effects$`scale(mean_ODBA)` %>% data.frame
m2_tenure_eff = m2_effects$`scale(log(tenure_days))` %>% data.frameFig: ODBA and mating success
success_odba_plot <- ggplot(ind_fertile_activity)+
geom_line(data = m2_active_eff,
aes(y = fit, x = mean_ODBA),
colour=col_scale[1]) +
geom_ribbon(data = m2_active_eff,
aes(y = fit, x = mean_ODBA, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[1])+
geom_jitter(aes(x=mean_ODBA, y=n_nests, shape=likely_extra_nests),
fill=male_colour,
colour=male_colour,
size=point_size,
alpha=point_alpha,
width=0.005, height=0.04) +
theme_classic()+
scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
scale_x_continuous(expand = c(0, 0), breaks=c(seq(0.1,0.5,by=0.1))) +
scale_y_continuous(breaks=c(seq(0,2,by=1))) +
coord_cartesian(xlim=c(0.1,0.53), ylim=c(0,2.1))+
xlab("\n mean ODBA (g)") +
ylab("mating success\n")+
theme(legend.position="right",
plot.margin=unit(c(0,0,0,0), 'lines'),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
success_odba_plot
Reproductive success (number of offspring sired)
Prediction 5 Increased time spent active (%) is associated with increased reproductive success Prediction 6 Males with longer site tenure have higher reproductive success Prediction 7 Males with longer site tenure have higher reproductive success
Model: % active and reproductive success
success_m3 <- glmmTMB(n_eggs ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day),
data=ind_fertile_activity,
# data=ind_fertile_activity %>% filter(likely_extra_nests=="less likely"),
family=nbinom2)
summary(success_m3)## Family: nbinom2 ( log )
## Formula:
## n_eggs ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
##
## AIC BIC logLik deviance df.resid
## 158.8 170.7 -74.4 148.8 76
##
##
## Dispersion parameter for nbinom2 family (): 0.149
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.611 0.335 -1.82 0.068 .
## scale(prop_activity) 0.306 0.416 0.74 0.462
## scale(log(tenure_days)) 0.603 0.435 1.38 0.167
## scale(arrival_day) -0.026 0.301 -0.09 0.931
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1simulationOutput <- simulateResiduals(fittedModel = success_m3, plot=F)
plot(simulationOutput)
success_m3_flex <- as_flextable(success_m3)
if(save_tables==TRUE){
save_as_docx(success_m3_flex,
pr_section = sect_properties,
path = glue("{output_path}/activity_reproduction_model.docx"))
}Similar effects are observed when males near predated nests are excluded.
m3_effects = effects::allEffects(success_m3, xlevels = 100)
m3_active_eff = m3_effects$`scale(prop_activity)` %>% data.frame
m3_tenure_eff = m3_effects$`scale(log(tenure_days))` %>% data.framem3_effects = as.data.frame(confint(success_m3, level=0.95))
write.csv(m3_effects,
file=glue("{output_supp_effects}/activity-repr-effects_{window_length}min_{data_type}-{threshold_type}thresh.csv"))Fig: % active and reproductive success
egg_active_plot <- ggplot(ind_fertile_activity)+
geom_line(data = m3_active_eff,
aes(y = fit, x = prop_activity*100),
colour=col_scale[2]) +
geom_ribbon(data = m3_active_eff,
aes(y = fit, x =prop_activity*100, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[2])+
geom_jitter(aes(x=prop_activity*100, y=n_eggs, shape=likely_extra_nests),
fill=male_colour,
colour=male_colour,
size=point_size,
alpha=point_alpha,
width=0.005, height=0.15) +
theme_classic()+
scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,90,by=20))) +
coord_cartesian(xlim=c(10,90), ylim=c(0,6.2))+
xlab("\n activity (% time)") +
ylab("reproductive success\n")+
theme(legend.position="right",
plot.margin=unit(c(0,0,0,0), 'lines'),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
egg_active_plot
Model: Tenure and reproductive success
success_mt2 <- glmmTMB(n_eggs ~ scale(log(tenure_days_exact)) + scale(arrival_day),
data=ind_tenure,
# data=ind_tenure %>% filter(potential_n_extra_nests==0),
family=nbinom2)
summary(success_mt2)## Family: nbinom2 ( log )
## Formula: n_eggs ~ scale(log(tenure_days_exact)) + scale(arrival_day)
## Data: ind_tenure
##
## AIC BIC logLik deviance df.resid
## 158.7 169.1 -75.3 150.7 95
##
##
## Dispersion parameter for nbinom2 family (): 0.139
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3452 0.4366 -3.08 0.0021 **
## scale(log(tenure_days_exact)) 1.9216 0.7530 2.55 0.0107 *
## scale(arrival_day) -0.0294 0.3356 -0.09 0.9301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(success_mt2, level=0.95)## 2.5 % 97.5 % Estimate
## (Intercept) -2.201 -0.489 -1.3452
## scale(log(tenure_days_exact)) 0.446 3.397 1.9216
## scale(arrival_day) -0.687 0.628 -0.0294simulationOutput <- simulateResiduals(fittedModel = success_mt2, plot=F)
plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.3, p-value = 0.4
## alternative hypothesis: two.sidedSame conclusions when males near predated nests are excluded.
success_mt2_flex <- as_flextable(success_mt2)
if(save_tables==TRUE){
save_as_docx(success_mt2_flex,
pr_section = sect_properties,
path = glue("{output_path}/tenure_reproduction_model.docx"))
}mt2_effects = effects::allEffects(success_mt2, xlevels = 100)
mt2_tenure_eff = mt2_effects$`scale(log(tenure_days_exact))` %>% data.frameFig: Tenure and reproductive success
egg_tenure_plot <- ggplot(ind_tenure) +
geom_line(data = mt2_tenure_eff,
aes(y = fit, x = tenure_days_exact),
colour=col_scale[1]) +
geom_ribbon(data = mt2_tenure_eff,
aes(y = fit, x = tenure_days_exact, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[1])+
geom_jitter(aes(x=tenure_days_exact, y=n_eggs, shape=likely_extra_nests),
fill=male_colour,
colour=male_colour,
size=point_size,
alpha=point_alpha,
width=0.3, height=0.15) +
theme_classic()+
scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,30,by=10))) +
coord_cartesian(xlim=c(0,32), ylim=c(0,6.2))+
xlab("\n tenure (days)") +
ylab("reproductive success\n") +
theme(legend.position="right",
plot.margin=unit(c(0,0,0,0), 'lines'),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
egg_tenure_plot
Model: ODBA and reproductive success
success_m4 <- glmmTMB(n_eggs ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day),
data=ind_fertile_activity,
# data=ind_fertile_activity %>% filter(likely_extra_nests=="less likely"),
family=nbinom2)
summary(success_m4)## Family: nbinom2 ( log )
## Formula:
## n_eggs ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
##
## AIC BIC logLik deviance df.resid
## 157.2 169.2 -73.6 147.2 76
##
##
## Dispersion parameter for nbinom2 family (): 0.163
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.692 0.331 -2.09 0.037 *
## scale(mean_ODBA) 0.629 0.425 1.48 0.139
## scale(log(tenure_days)) 0.483 0.418 1.16 0.247
## scale(arrival_day) 0.127 0.329 0.39 0.699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1simulationOutput <- simulateResiduals(fittedModel = success_m3, plot=F)
plot(simulationOutput)
Similar effects when males near predated nests excluded.
success_m4_flex <- as_flextable(success_m4)
if(save_tables==TRUE){
save_as_docx(success_m4_flex,
pr_section = sect_properties,
path = glue("{output_path}/odba_reproduction_model.docx"))
}m4_effects = effects::allEffects(success_m4, xlevels = 100)
m4_active_eff = m4_effects$`scale(mean_ODBA)` %>% data.frame
m4_tenure_eff = m4_effects$`scale(log(tenure_days))` %>% data.frameFig: ODBA and reproductive success
egg_odba_plot <- ggplot(ind_fertile_activity)+
geom_line(data = m4_active_eff,
aes(y = fit, x = mean_ODBA),
colour=col_scale[2]) +
geom_ribbon(data = m4_active_eff,
aes(y = fit, x = mean_ODBA, ymin = lower, ymax = upper),
alpha = 0.2, fill=col_scale[2])+
geom_jitter(aes(x=mean_ODBA, y=n_eggs, shape=likely_extra_nests),
fill=male_colour,
colour=male_colour,
size=point_size,
alpha=point_alpha,
width=0.005, height=0.15) +
theme_classic()+
scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,0.5,by=0.1))) +
coord_cartesian(xlim=c(0.1,0.52), ylim=c(0,6.2))+
xlab("\n mean ODBA (g)") +
ylab("reproductive success\n")+
theme(legend.position="right",
plot.margin=unit(c(0,0,0,0), 'lines'),
text=element_text(family=font.family),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size))
egg_odba_plot
Combined figure: Mating success and behaviour
fig3 <- success_active_plot +
success_odba_plot + ylab("\n") +
success_tenure_plot + ylab("\n") +
plot_layout(nrow = 1,
widths=c(1,1.05,1.05),
heights=c(1))+
plot_annotation(tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')') &
theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
plot.tag.position = c(.3, .96))
fig3
if(save_plots==TRUE){
ggsave(plot=fig3,
filename='figure3-beh-mating-success.jpg',
path=output_path,
width=16,
height=7,
units="cm",
dpi=300)
}Combined supplementary figure: Reproductive success and behaviour
supfigsuccess <- egg_active_plot +
egg_odba_plot + ylab("\n") +
egg_tenure_plot + ylab("\n") +
plot_layout(nrow = 1,
widths=c(1,1.05,1.05),
heights=c(1))+
plot_annotation(tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')') &
theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
plot.tag.position = c(.3, .96))
supfigsuccess
if(save_plots==TRUE){
ggsave(plot=supfigsuccess,
filename='figure-rsuccess.jpg',
path=output_path,
width=16,
height=7,
units="cm",
dpi=300)
}rm(m1_active_eff, m1_active_eff2, m1_effects, m1_effects2, m1_tenure_eff,
m2_active_eff, m2_effects, m2_tenure_eff,
m3_active_eff, m3_effects, m3_tenure_eff,
m4_active_eff, m4_effects, m4_tenure_eff,
success_m1, success_m1_flex,
success_m2, success_m2_flex,
success_m3, success_m3_flex,
success_m4, success_m4_flex,
success_mt, success_mt_flex)Q: Does activity relative to rivals predict success at a given nest?
For this analysis, we don’t want to limit the dataset just to males with at least 2 days of data.
But we do want to exclude nests that were depredated before genotyping.
We adjust the dataset accordingly.
Median daily distance from the nest site tended to be within 200 metres (in 14/18 cases).
The table below provides a summary of how many days each male who sired offspring was recorded, during the days when the female he mated with was fertile.
nest_male_summary <- nest_male_date %>%
mutate(dist_m=as.numeric(dist_m),
success=as.integer(success)) %>%
filter(success==1) %>%
select(nest, date, bird_ID, success, colour_combo) %>%
group_by(nest, bird_ID, success, colour_combo) %>%
summarise(n_days_father_recorded=n()) %>%
arrange(n_days_father_recorded) %>%
as.data.frame()
nest_male_summary <- nest_dates %>%
group_by(nest) %>%
summarise(n_days_mother_fertile=n()) %>%
full_join(., nest_parentage, by='nest') %>%
filter(father1!="unbanded") %>%
full_join(., nest_male_summary, by='nest') %>%
select(nest, colour_combo, n_days_mother_fertile, n_days_father_recorded) %>%
mutate(father_unrecorded_days=n_days_mother_fertile-n_days_father_recorded) %>%
arrange(father_unrecorded_days)
nest_male_summary %>% flextable()nest | colour_combo | n_days_mother_fertile | n_days_father_recorded | father_unrecorded_days |
|---|---|---|---|---|
P1401 | Y,Y,Y | 8 | 8 | 0 |
P203 | DB,R,O | 7 | 7 | 0 |
P214 | O,DG,DG | 8 | 8 | 0 |
P1003 | O,DB,O | 8 | 7 | 1 |
P208 | O,Y,O | 8 | 7 | 1 |
P401 | PI,PI,R | 8 | 7 | 1 |
P403 | PI,R,PI | 8 | 7 | 1 |
P1905 | DG,O,O | 8 | 6 | 2 |
P201 | Y,DG,PI | 8 | 5 | 3 |
P1903 | DB,DG,Y | 8 | 4 | 4 |
P701 | PI,PI,R | 8 | 4 | 4 |
P204 | DB,R,R | 8 | 3 | 5 |
P403 | O,DG,DG | 8 | 3 | 5 |
P702 | Y,Y,R | 8 | 2 | 6 |
P1001 | DG,DB,DG | 8 | 1 | 7 |
P1103 | DG,R,R | 8 | 1 | 7 |
P1403 | Y,DG,PI | 8 | 1 | 7 |
P701 | R,DG,PI | 8 | 1 | 7 |
nest_male_summary %>%
summarise(mean_unrecorded_days=mean(father_unrecorded_days),
sd_unrecorded_days=sd(father_unrecorded_days))## # A tibble: 1 × 2
## mean_unrecorded_days sd_unrecorded_days
## <dbl> <dbl>
## 1 3.39 2.70Model: Activity relative to competitors
nest_male_date_sub <- nest_male_date %>%
filter(!nest %in% predated_nests$nest) %>%
mutate(dist_m=as.numeric(dist_m),
success=as.integer(success)) %>%
select(nest, date, bird_ID, colour_combo, mean_ODBA, prop_activity, dist_m, success) %>%
filter(dist_m <= 200) %>%
group_by(date, nest) %>%
mutate(z_prop_act = scale(prop_activity),
z_ODBA = scale(mean_ODBA)) %>%
filter(!is.na(z_prop_act)) %>%
as.data.frame()nest_date_n <- nest_male_date_sub %>%
group_by(nest, date) %>%
summarise(n=n()) %>%
arrange(n)
hist(nest_date_n$n, breaks=20)
Median sample size (number of rival males per nest per day):
median(nest_date_n$n)## [1] 4Mean sample size:
mean(nest_date_n$n)## [1] 4.48Standard deviation:
sd(nest_date_n$n)## [1] 2.14plot(nest_male_date_sub$dist_m, nest_male_date_sub$z_prop_act)
nest_male_m <- glmmTMB(success ~ scale(z_prop_act) +
(1|bird_ID) +
(1|nest),
data = nest_male_date_sub,
family = binomial())
summary(nest_male_m)## Family: binomial ( logit )
## Formula: success ~ scale(z_prop_act) + (1 | bird_ID) + (1 | nest)
## Data: nest_male_date_sub
##
## AIC BIC logLik deviance df.resid
## 254 273 -123 246 987
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## bird_ID (Intercept) 117 10.8
## nest (Intercept) 114 10.7
## Number of obs: 991, groups: bird_ID, 62; nest, 40
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.179 4.671 -4.32 0.000016 ***
## scale(z_prop_act) 0.734 0.311 2.36 0.018 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1simulationOutput <- simulateResiduals(fittedModel = nest_male_m, plot=F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1, p-value = 0.3
## alternative hypothesis: two.sidedtestDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.6, p-value = 0.2
## alternative hypothesis: two.sidednm_rel_effects <- confint(nest_male_m)
nm_rel_effects <- cbind(name=rownames(nm_rel_effects),
data.frame(nm_rel_effects,
row.names=NULL)) %>%
rename(`2.5`=`X2.5..`,
`97.5`=`X97.5..`,
estimate=Estimate) %>%
mutate(effect_type="rel")if(save_tables==TRUE){
nest_m_flex <- as_flextable(nest_male_m)
save_as_docx(nest_m_flex,
pr_section = sect_properties,
path = glue("{output_path}/scaled_nest_male_model.docx"))
}Fig: Activity relative to competitors
nest_male_plotdat <- rbind(nest_male_date_sub %>% mutate(plot_type="box"),
nest_male_date_sub %>% mutate(plot_type="point"))
nest_male_plotdat$plot_factor <- factor(
paste0("s", nest_male_plotdat$success, "-", nest_male_plotdat$plot_type),
levels = c("s0-point", "s0-box", "s1-box", "s1-point"))
nest_male_boxplot <- ggplot() +
theme_classic()+
geom_point(data = nest_male_plotdat %>%
filter(plot_factor == "s0-point"),
aes(x = plot_factor,
y = z_prop_act),
size=small.point.size,
alpha=point_alpha,
shape=default_shape,
fill=male_colour,
position=position_jitter(w=0.08, h=0, seed = 1),
) +
geom_boxplot(data = nest_male_plotdat %>%
filter(plot_factor == "s0-box"),
aes(x = plot_factor,
y = z_prop_act),
outlier.shape=default_shape,
fill="white",
width=0.7) +
geom_boxplot(data = nest_male_plotdat %>%
filter(plot_factor == "s1-box"),
aes(x = plot_factor,
y = z_prop_act),
outlier.shape=default_shape,
fill="white",
width=0.7)+
geom_point(data = nest_male_plotdat %>%
filter(plot_factor == "s1-point"),
aes(x = plot_factor,
y = z_prop_act),
size=small.point.size,
alpha=point_alpha,
shape=default_shape,
fill=male_colour,
position=position_jitter(w=0.08, h=0, seed = 1),
) +
coord_cartesian(ylim=c(-2.8,2.8))+
scale_y_continuous(breaks=c(seq(-2,2,by=1)),
expand=c(0,0))+
scale_x_discrete(breaks = c("s0-point", "s1-point"),
labels = c("unsuccessful", "successful"))+
xlab("\n mating success")+
ylab("male activity\n (scaled relative to local males)\n")+
theme(text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
legend.position="right",
legend.justification="top",
legend.background = element_blank(),
legend.box.background = element_rect(colour = "grey"))
nest_male_boxplot
Model: ODBA relative to competitors
nest_male_m2 <- glmmTMB(success ~ scale(z_ODBA) +
(1|bird_ID) +
(1|nest),
data = nest_male_date_sub,
family = binomial())
summary(nest_male_m2)## Family: binomial ( logit )
## Formula: success ~ scale(z_ODBA) + (1 | bird_ID) + (1 | nest)
## Data: nest_male_date_sub
##
## AIC BIC logLik deviance df.resid
## 252 272 -122 244 987
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## bird_ID (Intercept) 120 11
## nest (Intercept) 121 11
## Number of obs: 991, groups: bird_ID, 62; nest, 40
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.575 4.577 -4.50 0.0000069 ***
## scale(z_ODBA) 0.912 0.353 2.58 0.0098 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1simulationOutput <- simulateResiduals(fittedModel = nest_male_m2, plot=F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1, p-value = 0.3
## alternative hypothesis: two.sidedtestDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.6, p-value = 0.3
## alternative hypothesis: two.sidednm_rel_effects2 <- confint(nest_male_m2)
nm_rel_effects2 <- cbind(name=rownames(nm_rel_effects2),
data.frame(nm_rel_effects2,
row.names=NULL)) %>%
rename(`2.5`=`X2.5..`,
`97.5`=`X97.5..`,
estimate=Estimate) %>%
mutate(effect_type="rel")if(save_tables==TRUE){
nest_m_flex <- as_flextable(nest_male_m2)
save_as_docx(nest_m_flex,
pr_section = sect_properties,
path = glue("{output_path}/scaled_nest_male_odba_model.docx"))
}Fig: ODBA relative to competitors
nest_male_boxplot2 <- ggplot() +
theme_classic()+
geom_point(data = nest_male_plotdat %>%
filter(plot_factor == "s0-point"),
aes(x = plot_factor,
y = z_ODBA),
size=small.point.size,
alpha=point_alpha,
shape=default_shape,
fill=male_colour,
position=position_jitter(w=0.08, h=0, seed = 1),
) +
geom_boxplot(data = nest_male_plotdat %>%
filter(plot_factor == "s0-box"),
aes(x = plot_factor,
y = z_ODBA),
outlier.shape=default_shape,
fill="white",
width=0.7) +
geom_boxplot(data = nest_male_plotdat %>%
filter(plot_factor == "s1-box"),
aes(x = plot_factor,
y = z_ODBA),
outlier.shape=default_shape,
fill="white",
width=0.7)+
geom_point(data = nest_male_plotdat %>%
filter(plot_factor == "s1-point"),
aes(x = plot_factor,
y = z_ODBA),
size=small.point.size,
alpha=point_alpha,
shape=default_shape,
fill=male_colour,
position=position_jitter(w=0.08, h=0, seed = 1),
) +
coord_cartesian(ylim=c(-2.8,2.8))+
scale_y_continuous(breaks=c(seq(-2,2,by=1)),
expand=c(0,0))+
scale_x_discrete(breaks = c("s0-point", "s1-point"),
labels = c("unsuccessful", "successful"))+
xlab("\n mating success")+
ylab("male ODBA\n (scaled relative to local males)\n")+
theme(text=element_text(family=font.family),
axis.title = element_text(size=axis.heading.size),
axis.text = element_text(size=label.size),
legend.title = element_text(size=legend.heading.size),
legend.text = element_text(size=label.size),
legend.position="right",
legend.justification="top",
legend.background = element_blank(),
legend.box.background = element_rect(colour = "grey"))
nest_male_boxplot2
Combined figure: Activity and ODBA relative to competitors
fig4 <- nest_male_boxplot +
plot_spacer() +
nest_male_boxplot2 +
plot_layout(nrow = 1,
widths=c(1,0.1,1),
heights=c(1))+
plot_annotation(tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')') &
theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
plot.tag.position = c(.28, .96))
fig4
if(save_plots==TRUE){
ggsave(plot=fig4,
filename='figure4-relative-male-activity.jpg',
path=output_path,
width=16,
height=7,
units="cm",
dpi=300)
}Q: How repeatable is activity during the fertile period?
repeat_lme <- lmer(prop_activity ~ (1|bird_ID),
data = daily_activity %>% filter(breeding_stage=="fertile"))
summary(repeat_lme)## Linear mixed model fit by REML ['lmerMod']
## Formula: prop_activity ~ (1 | bird_ID)
## Data: daily_activity %>% filter(breeding_stage == "fertile")
##
## REML criterion at convergence: -313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7510 -0.5959 -0.0531 0.5941 2.5074
##
## Random effects:
## Groups Name Variance Std.Dev.
## bird_ID (Intercept) 0.0169 0.130
## Residual 0.0260 0.161
## Number of obs: 555, groups: bird_ID, 81
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.6055 0.0162 37.4simulationOutput <- simulateResiduals(fittedModel = repeat_lme, plot=F)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 0.9
## alternative hypothesis: two.sided# Get variance
vc = VarCorr(repeat_lme)
withingroup_var=attr(vc, 'sc')[1]^2
betweengroup_var=attr(vc$bird_ID, 'stddev')[1]^2
# repeatability
R = betweengroup_var/(betweengroup_var+withingroup_var)
R## (Intercept)
## 0.394Within-individual repeatability of daily activity during the fertile phase was R = 0.394.